\documentclass[11pt]{article}       
\usepackage{a4wide,amsmath}
\usepackage{graphicx}

\setlength{\textheight}{25cm}
\setlength{\textwidth}{15cm}
\oddsidemargin=0mm
\evensidemargin=0mm
%\setlength{\parindent}{0mm}
\setlength{\itemsep}{-5pt plus-5pt minus-5pt}
\setlength{\parskip}{4pt plus1pt minus1pt}
\setlength{\parsep}{2pt plus1pt minus1pt}
\voffset -1.5cm

    \def\abl#1#2{\frac {d #1}{d #2}}
    \def\pabl#1#2{\frac {\partial #1}{\partial #2}}
    \def\eg{e.\,g.\ }
    \def\ie{i.\,e.\ }

\begin{document}

\sloppy

\centerline{\huge\bf Lambda Operator Technique in ProDiMo}
{\ }\\[-1.7ex]
\centerline{\large Peter Woitke, Feb.-March 2011}

\section{Source Function and Opacity Interpolation}

Given a point $\vec{x}\!=\!\{x,y,z\}$ in 3D, the source function
$S_\nu$ and extinction coefficient $\kappa_\nu^{\rm ext}$ are
linearly interpolated from the four nearest points in the grid in cylinder
coordinates $r\!=\!\sqrt{x^2+y^2}$ and $z'\!=\!|z|$
\begin{eqnarray}
  S_\nu(\vec{x}) &=& f_{i,j}   S_\nu^{i,j} 
                 + f_{i,j+1}   S_\nu^{i,j+1} 
                 + f_{i+1,j}   S_\nu^{i+1,j} 
                 + f_{i+1,j+1} S_\nu^{i+1,j+1}  \label{Snu}\\
  \kappa_\nu^{\rm ext}(\vec{x}) &=& f_{i,j}     \kappa^{\rm ext}_{i,j} 
                 + f_{i,j+1}   \kappa^{\rm ext}_{i,j+1} 
                 + f_{i+1,j}   \kappa^{\rm ext}_{i+1,j} 
                 + f_{i+1,j+1} \kappa^{\rm ext}_{i+1,j+1}
\end{eqnarray}
where $S_\nu^{i,j}$ and $\kappa^{\rm ext}_{i,j}$ are the pre-calculated
source functions and opacities at $\{r_i,z'_j\}$, see {\tt WHEREAMI()}
and {\tt INTERPOL\_RT()}.

\section{Ray integration}

Given a ray direction $\vec{n}$, the radiative transfer equation from
point $\vec{x}\!=\!\vec{x_0}+s\,\vec{n}$ to point $\vec{x_0}+(s+\Delta
s)\,\vec{n}=\vec{x}+\Delta s\,\vec{n}$ is to be solved based on two
points
\begin{eqnarray}
  S_1 &=& S_\nu(\vec{x}) \\
  S_2 &=& S_\nu(\vec{x}+\Delta s\,\vec{n}) \\
  \kappa_1 &=& \kappa_\nu^{\rm ext}(\vec{x}) \\
  \kappa_2 &=& \kappa_\nu^{\rm ext}(\vec{x}+\Delta s\,\vec{n}) 
\end{eqnarray}
The formal solution of the radiative transfer equation at point $\vec{x}_0$
backward along the ray $\vec{n}$ is
\begin{eqnarray}
  \tau_\nu(s) &=& \int_0^s \kappa_\nu^{\rm ext}(s')\, ds' \\
  I_\nu(\vec{x}_0,-\vec{n}) &=& 
     {\rm e}^{-\tau_\nu(s_{\rm max})}\,I_\nu^{\rm inc}
     + \int\limits_{s=0}^{s_{\rm max}} 
       {\rm e}^{-\tau_\nu(s)}\,S_\nu(s)\,\kappa_\nu^{\rm ext}(s)\,ds
  \label{formSol}
\end{eqnarray}
To arrive at a simple analytic expression for the solution of
Eq.(\ref{formSol}) we assume that $\kappa_\nu^{\rm ext}$ and $S_\nu$
are linear in $s$ for a small step $\Delta s$:
\begin{eqnarray}
  \kappa_\nu^{\rm ext}(s) &\approx& A+B\,s\\
  S_\nu(s) &\approx& C+D\,s
\end{eqnarray}
with
\begin{eqnarray}
  A &=& \kappa_1\\
  B &=& (\kappa_2-\kappa_1)/\Delta s\\
  C &=& S_1 \label{C}\\
  D &=& (S_2-S_1)/\Delta s \label{D}
\end{eqnarray}
The increments for a step $s \to s+\Delta s$ are
\begin{eqnarray}
  \tau_\nu(s+\Delta s) &=& \tau_\nu(s) 
                  + \int\limits_0^{\Delta s} (A+B\,s')\,ds'\nonumber\\
                  &=& \tau_\nu(s) + A\Delta s + \frac{B}{2}\Delta s^2 \\
                  &=& \tau_\nu(s) + \frac{1}{2}(\kappa_1+\kappa_2)\Delta s\\
  I_\nu(s+\Delta s) &=& I_\nu(s) + {\rm e}^{-\tau_\nu(s)}
          \int\limits_0^{\Delta s}
          {\rm e}^{-As' - \frac{B}{2}s'^2} (C+D\,s') (A+B\,s')\,ds'
  \label{formSol1}
\end{eqnarray}
The analytic solution of Eq.~(\ref{formSol1}) is quite awful and slow,
involving the error function. It is probably more efficient to make
more, smaller steps with further simplifications. We assume that
$\kappa_\nu^{\rm ext}(s)$ is constant with $\kappa_\nu^{\rm
ext}(s)\approx\bar{\kappa}=\frac{1}{2}(\kappa_1+\kappa_2)$, in
which case
\begin{eqnarray}
  \tau_\nu(s+\Delta s) &=& \tau_\nu(s)+ \bar{\kappa}\,\Delta s\\
  I_\nu(s+\Delta s) &=& I_\nu(s) + {\rm e}^{-\tau_\nu(s)}
          \int\limits_0^{\Delta s}
          \exp(-\bar{\kappa}s') (C+D\,s') \,\bar{\kappa}\,ds' \\
     &=& I_\nu(s) + {\rm e}^{-\tau_\nu(s)} \int\limits_0^t \exp(-x)
         \left(S_1+(S_2-S_1)\frac{x}{t}\right) \,dx \\
     &=& I_\nu(s) + {\rm e}^{-\tau_\nu(s)}
         \left[S_1+\frac{S_2-S_1}{t} 
         - \exp(-t)
         \left(S_1+\frac{S_2-S_1}{t}(1+t)\right)\right]\\
     &=& I_\nu(s) + {\rm e}^{-\tau_\nu(s)}
         \left[S_1\left(1 - \frac{1-{\rm e}^{-t}}{t}\right)
              +S_2\left(-{\rm e}^{-t}
                          + \frac{1-{\rm e}^{-t}}{t}\right)\right]
  \label{formSol2}
\end{eqnarray}
where $t=\bar{\kappa}\,\Delta s$. Equations~(\ref{C}) and Eq.~(\ref{D})
have been used.

\subsection{Numerical diffusion}

It is interesting to study the limiting behavior of the weights for $S_1$ and 
$S_2$ in the case $t\to\infty$, i.e. for a large single step in optical depth 
\begin{equation}
  I_\nu(s+\Delta s) \to I_\nu(s) + {\rm e}^{-\tau_\nu(s)}
                  \left[S_1\left(1 - \frac{1}{t}\right)
                       +S_2\left(\frac{1}{t}\right)\right] \ .
\end{equation}
The influence of neighboring source functions hence scales with
$1/\tau$, which is bad, a large influence, a problem if there are
large gradients in $S(\tau)$. To demonstrate this effect, let's assume
that $S(t)\propto t^2$ between $S_1$ and $S_2$
\begin{eqnarray}
 && I_\nu(s+\Delta s) = I_\nu(s) + {\rm e}^{-\tau_\nu(s)}
          \int\limits_0^t \exp(-x)
          \left(S_1+(S_2-S_1)\left(\frac{x}{t}\right)^2\right) \,dx 
  \label{SquareDep}\\
  &=& I_\nu(s) + {\rm e}^{-\tau_\nu(s)} 
         \Bigg[S_1\left(1-\frac{2}{t^2}\Big(1-{\rm e}^{-t}(1+t)\Big)\right)
 +S_2\left(-{\rm e}^{-t} +\frac{2}{t^2}\Big(1-{\rm e}^{-t}(1+t)\Big)
       \right)\Bigg] \nonumber
\end{eqnarray}
In this case, the influence of neighboring source functions at large
optical depth would scale with $1/\tau^2$ (much less, i.e.\ much
better). Unfortunately, the approach (\ref{SquareDep}) is inconsistent
with the way that source functions are interpolated within the grid,
which uses a {\it linear} interpolation between the 4 nearest points
in the grid (see Eq.~\ref{Snu}).

The large numerical diffusion was the reason why ProDiMo up to revision 
995 used a logarithmic interpolation of the source function in the first place,
which is, however, inconsistent with the idea of a linear $\Lambda$ operator.
Thus, I presently do not have a good idea how to build a precise and fast
converging radiative transfer method for ProDiMo. The old scheme is more
precise, but does not converge for large optical depth -- the new schema
does converge but has large numerical diffusion.


\section{The $\Lambda$ operator}

In the following, the idea is to express $I_\nu$ by a linear superposition of
source functions as
\begin{equation}
  I_\nu(\vec{x_0},-\vec{n}^{k,l}) = \alpha^{k,l} I_\nu^{\rm inc}
                   + \sum\limits_{ij} \lambda_{i,j}^{k,l} S_\nu^{i,j}
  \label{Ilambda}
\end{equation} 
which is possible because Eq.~(\ref{formSol2}) is {\it linear} in
$S_\nu$.  In contrast, in ProDiMo up to revision 995, a
log-interpolation is used for the source function. In that case, it is
not possible to derive Eq.~(\ref{Ilambda}).  $S_\nu^{i,j}$ refers to all
source functions at all grid points $\{r_i,z'_j\}$, although in
practise only a very few (the ones close to the ray $kl$) will
actually contribute to $ I_\nu(\vec{x_0},-\vec{n}_{kl})$.

For every step $\Delta s$, there are 8 entries $\lambda_{i,j}^{k,l}$ to
be incremented, 4 around the starting point $\vec{x}$ (with indices and
weights as $f_{i1,j1}$ ... $f_{i1+1,j1+1}$) and 4 around the end point
$\vec{x}+\Delta s\,\vec{n}^{kl}$ (with indices and weights as $f_{i2,j2}$
... $f_{i2+1,j2+1}$). In order to get these entries we have to use
Eq.~(\ref{Snu}) twice, and to combine it with Eq.~(\ref{formSol2}).
\begin{eqnarray}
\lambda_{i1,j1}^{k,l}     &\leftarrow& {\rm e}^{-\tau_\nu(s)} 
           \Big(1-\frac{1-{\rm e}^{-t}}{t}\Big) f_{i1,j1}     \nonumber\\ 
\lambda_{i1,j1+1}^{k,l}   &\leftarrow& {\rm e}^{-\tau_\nu(s)} 
           \Big(1-\frac{1-{\rm e}^{-t}}{t}\Big) f_{i1,j1+1}   \nonumber\\ 
\lambda_{i1+1,j1}^{k,l}   &\leftarrow& {\rm e}^{-\tau_\nu(s)} 
           \Big(1-\frac{1-{\rm e}^{-t}}{t}\Big) f_{i1+1,j1}   \nonumber\\ 
\lambda_{i1+1,j1+1}^{k,l} &\leftarrow& {\rm e}^{-\tau_\nu(s)} 
           \Big(1-\frac{1-{\rm e}^{-t}}{t}\Big) f_{i1+1,j1+1} \nonumber\\ 
\lambda_{i2,j2}^{k,l}     &\leftarrow& {\rm e}^{-\tau_\nu(s)}
 \Big(\frac{1}{t}-{\rm e}^{-t}(1+\frac{1}{t})\Big) f_{i2,j2}    \nonumber\\
\lambda_{i2,j2+1}^{k,l}   &\leftarrow& {\rm e}^{-\tau_\nu(s)}
 \Big(\frac{1}{t}-{\rm e}^{-t}(1+\frac{1}{t})\Big) f_{i2,j2+1}  \nonumber\\
\lambda_{i2+1,j2}^{k,l}   &\leftarrow& {\rm e}^{-\tau_\nu(s)}
 \Big(\frac{1}{t}-{\rm e}^{-t}(1+\frac{1}{t})\Big) f_{i2+1,j2}  \nonumber\\
\lambda_{i2+1,j2+1}^{k,l} &\leftarrow& {\rm e}^{-\tau_\nu(s)}
 \Big(\frac{1}{t}-{\rm e}^{-t}(1+\frac{1}{t})\Big) f_{i2+1,j2+1}
  \label{LamCoeff}
\end{eqnarray}
Integration of Eq.~(\ref{Ilambda}) over the solid angle, involving all 
ray directions $\vec{n}^{kl}$ with integration weights $\Delta\Omega_{k,l}$ 
yields
\begin{equation}
\boxed{
 J_\nu(\vec{x_0}) = \Lambda_\nu^\star(\vec{x_0}) I_\nu^\star 
                  + \Lambda_\nu^{\rm ISM}(\vec{x_0}) I_\nu^{\rm ISM}
    + \sum\limits_{i,j} \Lambda_\nu^{i,j}(\vec{x_0}) S_\nu^{i,j}
}
  \label{Jlambda}
\end{equation}
where the full $\Lambda$ operator coefficients are given by
\begin{equation}
  \Lambda_\nu^{i,j} = \frac{1}{4\pi} 
                \sum\limits_{k,l}\lambda_{i,j}^{k,l} \,\Delta\Omega_{k,l}\ .
\end{equation}
Similar expressions can be found for the incident coefficients
$\Lambda_\nu^\star$ and $\Lambda_\nu^{\rm ISM}$.

\subsection{Some properties of the $\Lambda$ operator}

Equation~(\ref{Jlambda}) states that every mean intensity in the grid
$J_\nu$ results from a linear superposition of all grid source functions
$S_\nu^{i,j}$ and the incident (stellar and ISM) intensities. The
coefficients $\Lambda_\nu^{i,j}$, $\Lambda_\nu^\star$,
$\Lambda_\nu^{\rm ISM}$ are given purely by grid geometry and optical
depths.  This is particularly practical in case of dust radiative
transfer, where the opacities, and hence the optical depths, are
temperature-independent. Here, the dust opacity structure fully
determines the coefficients in Eq.~(\ref{Jlambda}).  Thus, in
principle, the ray tracing needs to be done only once. All subsequent
iterations to solve the scattering problem and the search for the
proper dust temperature structure can be done directly with
Eq.~(\ref{Jlambda}). Unfortunately, storing the full $\Lambda$
operator is usually impossible. In 2D, the $\Lambda$ operator has rank
5, with typical total size as
\begin{equation}
 \rm (no.\ of\ grid\ points)^2 \times (no.\ of\ frequency\ points)
 \approx (100 \times 100)^2 \times 40 = 4\times 10^9
 \nonumber
\end{equation}
For real$\star$8 precision, this would require about 30\,GB of memory.

\section{The approximate $\Lambda$ operator}

The idea of the approximate $\Lambda$ operator is to store only a very
few elements of the full $\Lambda$ operator, namely the elements that
express the local feedbacks between neighboring grid points.  We omit
lower index $\nu$ for now, consider ${ik}$ as destination index and
$jl$ as source index, and abbreviate $\Lambda_\nu^{\star}=A$ and
$\Lambda_\nu^{\rm ISM}=C$. Then the full $\Lambda$ operator equation
(\ref{Jlambda}) and the partial derivatives are
\begin{eqnarray}
  J_{ik} &=& A_{ik} I_\star + C_{ik} I_{\rm ISM}
             + \sum\limits_{jl} \Lambda_{ik}^{jl} S_{jl} \label{fullLam}\\ 
  \frac{\partial{J_{ik}}}{{\partial S_{jl}}} &=& \Lambda_{ik}^{jl}
\end{eqnarray}
Making the leap to the approximate operator $\bar\Lambda$ means to throw
away most of the $\Lambda_{ik}^{jl}$ coefficients, and to store only the
nearest neighbor feedbacks:
\begin{equation}
  \bar\Lambda_{ik}^{jl} = \left\{\begin{array}{rl} 
                          \Lambda_{ik}^{jl} & \mbox{$ik$ ``close'' to $jl$}\\
                                          0 & \mbox{otherwise}
                          \end{array}\right.
\end{equation}
which means that
\begin{eqnarray}
  J_{ik} &\neq& A_{ik} I_\star + C_{ik} I_{\rm ISM}
             + \sum\limits_{jl} \bar\Lambda_{ik}^{jl} S_{jl} \\ 
  \frac{\partial{J_{ik}}}{{\partial S_{jl}}} &\approx& 
                          \bar\Lambda_{ik}^{jl} \label{AppLamOp}
\end{eqnarray}
The trick of the approximate $\bar\Lambda$ operator technique is to
use it for iterative corrections based on the fully linearized
problem.  The advantage of the $\bar\Lambda$ operator technique (in
contrast to more simple $\Lambda$ iteration techniques) is its fast
convergence even at high optical depths. In an optically thick
situation, the diagonal elements in the $\Lambda$ operator (how much
the local source function determines the mean intensity) are close to
1. The important energy exchange between different cells concerns
mainly neighboring cells (expressed by the entries in the secondary
diagonals), whereas the feedback of more distant cells becomes
negligible. Therefore, in the optically thick case, it is physically
sufficient to use the approximate (e.g.\ tri-diagonal) $\Lambda$
operator. The optically thin case converges rapidly anyway. In 2D,
there are maximum 8 cells around a central cell, hence the memory
consumption is
\begin{equation}
 \rm (no.\ of\ grid\ points) \times 9 \times (no.\ of\ frequency\ points)
 \approx (100 \times 100) \times 9 \times 40 = 3.6\times 10^6
 \nonumber
\end{equation}
which requires 27\,MB memory. 



\section{Applications}

\subsection{The scattering problem}

The scattering problem is that the calculation of the mean intensities
via a ray-tracing solution of the radiative transfer equation requires
to know the source functions, which in return depend on the mean
intensities.  For LTE and isotropic scattering we have
\begin{equation}
  S_\nu = \frac{\kappa_\nu^{\rm abs} B_\nu(T) + \kappa_\nu^{\rm sca} J_\nu}
               {\kappa_\nu^{\rm ext}}
        = (1-a_\nu) B_\nu(T) + a_\nu J_\nu
  \label{source}
\end{equation}
Let's assume that the dust temperatures are fixed and known for this
example. The $\Lambda$ iteration technique would start with a guess of
all $J_\nu$, then use Eq.~(\ref{source}) to calculate the source
functions, then perform a ray tracing (Eq.~\ref{formSol}) to calculate
$J_\nu$, which are again used as improved guess for Eq.~(\ref{source}),
and so on. At large optical depths, and large albedos $a_\nu =
\kappa_\nu^{\rm sca}/\kappa_\nu^{\rm ext}$, the scattering problem is
known to be a hard one for the $\Lambda$ iteration, with very slow convergence.

With the full $\Lambda$ operator technique, it would be possible to
solve the scattering problem in one step. We just have to combine
Eq.~(\ref{source}) with Eq.~(\ref{Jlambda}) and solve it. This can be
done separately for every frequency point.
\begin{eqnarray}
  J_{ik} &=& A_{ik} I^\star + C_{ik} I^{\rm ISM}
             + \sum\limits_{jl} \Lambda_{ik}^{jl} S_{jl} \nonumber\\
         &=& A_{ik} I^\star + C_{ik} I^{\rm ISM}
             + \sum\limits_{jl} \Lambda_{ik}^{jl} 
               \Big[\left(1-a_{jl}\right) B(T_{jl}) + a_{jl} J_{jl}\Big]
\end{eqnarray}

\begin{equation}
  \sum\limits_{jl} \Big(\delta_{ik}^{jl} - a_{jl}\Lambda_{ik}^{jl}\Big)J_{jl}
         = A_{ik} I^\star + C_{ik} I^{\rm ISM}
         + \sum\limits_{jl} \Lambda_{ik}^{jl} 
               \left(1-a_{jl}\right) B(T_{jl})
  \label{solveScat}
\end{equation}
The solution of the scattering problem Eq.~(\ref{solveScat}) is a linear
equation system, huge in dimension though. 

\subsubsection{Method 1}
Switching to the approximate $\bar\Lambda$ treatment, we perform one
ray tracing to get $J_\nu$ from $S_\nu$, and then consider Eq.~(\ref{source})
as a function $F_{ik}$ to be nullified. According to the Newton-Raphson
method, we are searching for corrections $DF\cdot \delta S_{jl} = -F_{ik}$ 
to improve the solution
\begin{eqnarray}
  F_{ik} &=& (1-a_{ik}) B(T_{ik}) + a_{ik} J_{ik} - S_{ik} \label{scat1}\\
  \frac{\partial F_{ik}}{\partial S_{jl}}&=&
             a_{ik}\,\bar\Lambda_{ik}^{jl} - \delta_{il}^{jl} \label{scat2}
\end{eqnarray} 
Equations (\ref{scat1}) and (\ref{scat2}) give the recipes for
one Newton-Raphson iteration step to improve the scattering solution.
Once the corrections $\delta S_{ik}$ are applied, the ray tracing needs
to be repeated.

\subsubsection{Method 2}

Although there is nothing wrong, in principle, with the scattering solution
discussed in the previous section, there are a few disadvantages of the method
described before
\begin{itemize}
\item only use of current $J_{ik}$
\item corrections to $S_{ik}$ instead of $J_{ik}$, which can result in
      negative $J_{ik}$
\end{itemize}
The following method makes use of two variants of $J_{ik}$, namely
its current status at iteration m, denoted as $J^{(m)}_{ik}$, and
the results after performing one formal solution, denoted as $J^{FS}_{ik}$.
A direct $\Lambda$-iteration would make corrections to $J^{(m)}_{ik}$ 
as
\begin{equation}
  \mbox{$\Lambda$-iteration:\quad} \delta J_{ik}\;=\;J^{FS}_{ik}-J^{(m)}_{ik}
\end{equation}
For our $\Lambda$-operator method we seek corrections $\delta J_{ik}$
that, after calculating the source functions and performing a formal
solution, would result in zero changes
\begin{eqnarray}
  J^{(m)}_{ik} + \delta J_{ik} 
    &=& A_{ik} I^\star + C_{ik} I^{\rm ISM}
      + \sum\limits_{jl} \Lambda_{ik}^{jl} 
        \Big[\left(1-a_{jl}\right) B(T_{jl}) 
             +a_{jl} \big(J^{(m)}_{jl}+\delta J_{jl}\big)\Big]\nonumber\\
    &=& A_{ik} I^\star + C_{ik} I^{\rm ISM}
      + \sum\limits_{jl} \Lambda_{ik}^{jl} 
        \Big[\left(1-a_{jl}\right) B(T_{jl}) +a_{jl} J^{(m)}_{jl}\Big]
      + \sum\limits_{jl} a_{jl}\,\Lambda_{ik}^{jl}\,\delta J_{jl}\nonumber\\
    &=& J^{FS}_{ik} 
      + \sum\limits_{jl} a_{jl}\,\Lambda_{ik}^{jl}\,\delta J_{jl}\nonumber\\
    &\approx& J^{FS}_{ik} 
      + \sum\limits_{jl} a_{jl}\,\bar\Lambda_{ik}^{jl}\,\delta J_{jl}
\end{eqnarray}
which results in
\begin{equation}
  \mbox{$\Lambda$-operator:\quad}
  \sum\limits_{jl} \big(\delta_{jl}^{kl} 
      - a_{jl}\,\bar\Lambda_{ik}^{jl}\big)\,\delta J_{jl} 
    \;=\; J^{FS}_{ik}-J^{(m)}_{ik}
  \label{scat3}
\end{equation}
Equation (\ref{scat3}) demonstrates that the $\Lambda$-operator
corrections are a kind of ``booster'' for the $\Lambda$-iteration
corrections. If scattering is unimportant $(a_{jl}\to 0)$, or if the
medium is optically thin $(\Lambda_{ik}^{jl}\to 0)$, the two methods
make identical corrections. In contrast, if the albedo is large
$(a_{jl}\to 1)$ and the medium is optically thick
$(\Lambda_{ik}^{jl}\to \delta_{jl}^{kl})$, then the $\Lambda$-operator
corrections become much larger as compared to the $\Lambda$-iteration
corrections.


\subsection{Determining the dust temperature structure}

\subsubsection{Method 1}
The energy balance of dust grains is written as
\begin{equation}
  \frac{\Gamma_{ik}}{4\pi} 
  + \sum_n \Big(J_{ik}^n - B^n(T_{ik})\Big)\kappa^n_{ik}\,\Delta\nu_n = 0
  \label{temp}
\end{equation}
$\Gamma_{ik}$ is a non-radiative heating rate of the dust, $\kappa^n_{ik}$
is the absorption coefficient at frequency point $\nu_n$ and spatial 
destination point $ik$. $\Delta\nu_n$ is the frequency integration weight.
In order to find the dust temperature structure, we have to solve
Eqs.~(\ref{temp}), (\ref{source}) and (\ref{Jlambda}). 

Due to the non-linear behavior of the Planck function, we cannot expect
to solve this problem in one step, but we need an iterative solution method,
for example a Newton-Raphson iteration. We consider Eq.~(\ref{temp}) as the
function $F_{ik}$ to be nullified, and need to work out the partial
derivatives $\partial{F_{ik}}/\partial{T_{jl}}$
\begin{equation}
  F_{ik} = \frac{\Gamma_{ik}}{4\pi} + \sum_n 
     \Big(J_{ik}^n - B^n(T_{ik})\Big)\kappa^n_{ik}\,\Delta\nu_n 
  \label{NewRaph1}
\end{equation}
By means of Eq.~(\ref{AppLamOp}) and the abbreviation
$B^n_{jl}=B^n(T_{jl})$ we find
\begin{equation}
  \frac{\partial F_{ik}}{\partial T_{jl}}
  = \sum_n \Big(\bar\Lambda_{ik}^{njl} 
                \pabl{S^n_{jl}}{T_{jl}}
                -\delta_{ik}^{jl} \pabl{B^n_{jl}}{T_{jl}} 
           \Big)\kappa^n_{ik}\,\Delta\nu_n 
\end{equation}
and inserting Eq.~(\ref{source}) yields
\begin{equation}
  \frac{\partial S^n_{jl}}{\partial T_{jl}} = 
   (1-a^n_{jl})\pabl{B^n_{jl}}{T_{jl}}
   + a^n_{jl}\frac{\partial J^n_{jl}}{\partial T_{jl}} \ .
\end{equation}
Using Eq.~(\ref{AppLamOp}) once more for the diagonal derivative 
$\partial J^n_{jl}/\partial T_{jl}$ we find
\begin{equation}
  \frac{\partial S^n_{jl}}{\partial T_{jl}} = 
   (1-a^n_{jl})\pabl{B^n_{jl}}{T_{jl}}
   + a^n_{jl}\bar\Lambda_{jl}^{njl} \pabl{S^n_{jl}}{T_{jl}}
\end{equation}
\begin{equation}
  \frac{\partial S^n_{jl}}{\partial T_{jl}}
  \left(1-a^n_{jl}\bar\Lambda_{jl}^{njl}\right) = 
   (1-a^n_{jl})\pabl{B^n_{jl}}{T_{jl}}
  \label{dSdT}
\end{equation}
which results in  
\begin{equation}
  \frac{\partial F_{ik}}{\partial T_{jl}}
  = \sum_n \left(\frac{\bar\Lambda_{ik}^{njl}\big(1-a^n_{jl}\big)}
                      {1-a^n_{jl}\bar\Lambda_{jl}^{njl}}
                -\delta_{ik}^{jl}\right)
           \pabl{B^n_{jl}}{T_{jl}} 
           \kappa^n_{ik}\,\Delta\nu_n 
  \label{NewRaph2}
\end{equation}
Equations (\ref{NewRaph1}) and (\ref{NewRaph2}) give the recipes for
one Newton-Raphson iteration step to improve the dust temperature
structure as $T_{ik}^{\rm new}=T_{ik}^{\rm old}+\delta T_{ik}$.

To embed the iterative solution scheme in ProDiMo, we must calculate
the effect of temperature changes $\delta T_{ik}$ on the mean
intensities, with exactly the same assumptions as used for the
calculation of $\delta T_{ik}$. 
\begin{eqnarray}
  \delta J_{ik}^{n} 
  &=& \sum_{jl} \bar\Lambda_{ik}^{njl} \delta S_{jl}^{n}\\
  &=& \sum_{jl} \bar\Lambda_{ik}^{njl} 
      \frac{\partial S_{jl}^{n}}{\partial T_{jl}} \delta T_{jl}
\end{eqnarray}
Using Eq.~(\ref{dSdT}) we find
\begin{equation}
  \delta J_{ik}^{n} = \sum_{jl} 
   \frac{\bar\Lambda_{ik}^{njl}\big(1-a^n_{jl}\big)}
                      {1-a^n_{jl}\bar\Lambda_{jl}^{njl}} 
   \pabl{B^n_{jl}}{T_{jl}} 
                 \delta T_{jl}
\end{equation}
These equations express the impact of {\it local} temperature changes
on the source function, and their global effect on the mean
intensities.  In order to calculate the total linear {\it global}
feedback between source functions and mean intensities, we must
perform, in addition, a scattering solution as described in the previous
section.

\newpage
\subsubsection{Method 2}

Method 1 is not fully satisfactory, because of the inconsistent
temperature and scattering corrections. We are now seeking simultaneous
mean intensity corrections $\delta J_{ik}^n$ and temperature
corrections $\delta T_{ik}$ such that
\begin{equation}
  \frac{\Gamma_{ik}}{4\pi} + \sum_n \Big(
     J_{ik}^n + \delta J_{ik}^n 
   - B_{ik}^n - \pabl{B_{ik}^n}{T_{ik}}\delta T_{ik} \Big)
     \kappa^n_{ik}\,\Delta\nu_n = 0
   \label{T2_equil}
\end{equation}
\centerline{and}
\begin{eqnarray}
  J^n_{ik} + \delta J^n_{ik} 
    &=& A^n_{ik} I_n^\star + C^n_{ik} I_n^{\rm ISM}
      + \sum\limits_{jl} \Lambda_{ik}^{njl} \left(S^n_{jl} + \delta S^n_{jl}\right)
      \nonumber\\
    &\approx& J^{n,FS}_{ik} 
      + \sum\limits_{jl} \Lambda_{ik}^{njl} \delta S^n_{jl} \label{T2a_RTsol}\\
    &=& J^{n,FS}_{ik} 
      + \sum\limits_{jl} \bar\Lambda_{ik}^{njl} 
        \Big[\left(1-a^n_{jl}\right) \pabl{B_{jl}^n}{T_{jl}}\delta T_{jl}
             +a^n_{jl} \delta J^n_{jl} \Big]\ . 
   \label{T2_RTsol}
\end{eqnarray}
Unfortunately, I don't see how to properly eliminate either $\delta J^n_{ik}$ or
$\delta T_{ik}$ from Eqs.~(\ref{T2_equil}) and (\ref{T2_RTsol}) to
deduce a linear system of equations for one of these corrections. In
the following, we try to solve Eqs.~(\ref{T2_equil}) and
(\ref{T2a_RTsol}) in an approximate way; we start be replacing
$J^n_{ik} + \delta J^n_{ik}$ in Eq.~(\ref{T2_equil}) by the r.h.s.
of Eq.~(\ref{T2a_RTsol})
\begin{equation}
  F_{ik} = \frac{\Gamma_{ik}}{4\pi} + \sum_n \bigg(
     J^{n,FS}_{ik} 
      + \sum\limits_{jl} \bar\Lambda_{ik}^{njl} \delta S^n_{jl}
      - B_{ik}^n - \pabl{B_{ik}^n}{T_{ik}}\delta T_{ik} \bigg)
     \kappa^n_{ik}\,\Delta\nu_n 
\end{equation}
Now, how does the energy balance at $ik$ change as we apply
a small T-change at $jl$\,? 
\begin{equation}
  \pabl{F_{ik}}{T_{jl}} = \sum_n \left(
      \bar\Lambda_{ik}^{njl} \pabl{\delta S^n_{jl}}{T_{jl}} 
   - \delta_{ik}^{jl}\pabl{B_{jl}^n}{T_{jl}} \right)
     \kappa^n_{ik}\,\Delta\nu_n
\end{equation}
Using Eq.~(\ref{dSdT}) yields
\begin{equation}
  \pabl{F_{ik}}{T_{jl}} = \sum_n \left(
   \frac{\bar\Lambda_{ik}^{njl}(1-a^n_{jl})}{1-a^n_{jl}\bar\Lambda_{jl}^{njl}} 
   - \delta_{ik}^{jl}\right) \pabl{B_{jl}^n}{T_{jl}}
     \kappa^n_{ik}\,\Delta\nu_n
\end{equation}
This way, we perform a temperature correction $\delta T_{ik}$. Given
the renewed temperartures, we determine new the mean intensity
corrections $\delta J^n_{ik}$ from Eq.~(\ref{T2a_RTsol}). 
\begin{equation}
  \sum\limits_{jl} \left(\delta_{ik}^{jl} - a^n_{jl}
                        \bar\Lambda_{ik}^{njl}\right) \delta J_{jl}^n 
  \;=\; \left( J^{n,FS}_{ik} - J^n_{ik}\right)
      + \sum\limits_{jl} \bar\Lambda_{ik}^{njl} 
        \left(1-a^n_{jl}\right) \pabl{B_{jl}^n}{T_{jl}}\delta T_{jl}
\end{equation}
The temperature and mean intensity corrections are then performed
iteratively (at fixed $\Lambda^{njl}_{ik}, a^n_{jl}$ that is), until convergence 
is achieved, before the next formal solution is carried out.




\newpage

\begin{figure}
\begin{tabular}{cc}
  PDR benchmark & Disk benchmark\\
  \hspace*{-11mm}\includegraphics[width=85mm]{epsfig/slab_LOPlin_resu.eps} 
 &\hspace*{ -7mm}\includegraphics[width=85mm]{epsfig/RT2_Td_LOPlin.eps}
\end{tabular}
\hspace*{0mm}\parbox{15.5cm}{\caption{Disappointing results from the
    $\Lambda$-operator technique when using a bi-linear interpolation
    of the source functions. It converges rapidly, but the results are
    dominated by numerical diffusion. The scattering problem for a
    semi-infinite slab (l.h.s.) results in much too large mean
    intensities. The calculated temperature structure for the disk RT
    benchmark (Pinte et al.~2009) is entirely wrong, much too high in
    the central midplane, probably due to radial numerical diffusion.
    }\label{fig1}}
\end{figure}

\section{A different approach}
The results of the $\Lambda$-operator treatment as described so far
are entirely disappointing. The reasons lie in the numerical
diffusion, which dominates over the proper physical interactions in
all cases we are interested in. Figure~\ref{fig1} shows two
examples. One way out of this dilemma is to switch to a
log-interpolation scheme for opacities and, in particular, for the
source functions as
\begin{eqnarray}
  \ln S_1 &\!\!\!=\!\!\!& f^{(1)}_{i1,k1}  \ln S_{i1,k1}
            + f^{(1)}_{i1,k1+1}   \ln S_{i1,k1+1} 
            + f^{(1)}_{i1+1,k1}   \ln S_{i1+1,k1} 
            + f^{(1)}_{i1+1,k1+1} \ln S_{i1+1,k1+1} \nonumber\\
  \ln S_2 &\!\!\!=\!\!\!& f^{(2)}_{i2,k2}  \ln S_{i2,k2}
            + f^{(2)}_{i2,k2+1}   \ln S_{i2,k2+1} 
            + f^{(2)}_{i2+2,k2}   \ln S_{i2+1,k2} 
            + f^{(2)}_{i2+2,k2+1} \ln S_{i2+1,k2+1} \nonumber \ ,
\end{eqnarray}
which, as already mentioned in the first ProDiMo paper (Woitke et
al.~2009), suffers much less from numerical diffusion. The disadvantage
of this approach is, however, that the influence of the source functions on the
RT results are not anymore linear -- the $\Lambda$-operator equation
(Eq.~\ref{fullLam}) is not valid anymore.

However, this is actually not required and Eq.~(\ref{fullLam}) is
never used in all practical applications, only as notation 
for the act of performing a formal RT solution. The actual use of the
$\bar\Lambda$-operator, in all applications, always considers
``small'', i.e. linear, deviations of the source functions are their
global feedback on the RT results.
\begin{equation}
  \bar\Lambda_{ik}^{jl} = \pabl{J_{ik}}{S_{jl}}
\end{equation}
where now, in general, the $\Lambda$-operator coefficients
$\Lambda_{ik}^{jl}$ depend on the source functions! The following
approach is hence similar to solving non-linear equations by means of
a complete linearisation of the problem and then making corrections based
on the linear dependencies, just like a Newton-Raphson iteration.


Assuming that both the extinction coefficient $\kappa_{\rm
ext}=\kappa_1\exp(Bs)$ and the source function $S=S_1\exp(Ds)$ 
vary log-linear along one transfer integration step we have
\begin{eqnarray}
  \tau_\nu(s+\Delta s) &=& \tau_\nu(s) 
        + \int\limits_0^{\Delta s} \kappa_1\exp(B\,s')\,ds'\nonumber\\
    &=& \tau_\nu(s) + \frac{\kappa_1}{B} \big(\exp(B\Delta s)-1\big)
    \label{fulltau}\\
  I_\nu(s+\Delta s) &=& I_\nu(s) + 
     \int\limits_0^{\Delta s} \kappa^{\rm ext} S \exp(-\tau_\nu) \,ds'
    \label{fullint}\\
         &=& I_\nu(s) + {\rm e}^{-\tau_\nu(s)} \,\zeta S_1 \kappa_1
     \int\limits_0^{\Delta s} \exp(B\,s' + D\,s' - \kappa_1\,s')\,ds'\nonumber\\
    &=& I_\nu(s) + {\rm e}^{-\tau_\nu(s)} \,\zeta \kappa_1 \Delta s \, S_1
                  \,\frac{1-\exp(-t)}{t}
\end{eqnarray}
with $t\!=\!(\kappa_1-B-D)\Delta s$. $\zeta$
is a correction factor close to 1, the ratio between the proper
integral (using Eq.~\ref{fulltau} for the optical depth in
Eq.~\ref{fullint}) and the ``fast'' integral when assuming
$\tau_\nu(s)\!=\kappa_1\cdot s$. After some analytic work using
$D=(\ln S_2-\ln S_1)/\Delta s$, we find
\begin{eqnarray}
  \pabl{}{S_1}\Big(S_1\frac{1-\exp(-t)}{t}\Big) 
    &=& \frac{1}{t}\Big(1-\frac{1-\exp(-t)}{t}\Big) \;=\;g_1(t) \\ 
  \pabl{}{S_2}\Big(S_1\frac{1-\exp(-t)}{t}\Big)
    &=& \frac{S_1}{S_2}\,
        \frac{1}{t}\left(\frac{1}{t}-\exp(-t)\Big(1+\frac{1}{t}\Big)\right)
                                    \;=\;\frac{S_1}{S_2}\,g_2(t) 
\end{eqnarray} 
In order to construct the $\Lambda$-operator entries we have to calculate
how the calculated intensities change as we vary the source functions
at the corner points of the log-interpolation:
\begin{eqnarray}
  \pabl{\Delta I_\nu}{S_{jl}} 
    &=& \pabl{\Delta I_\nu}{S_1} \pabl{S_1}{\ln S_1}
        \pabl{\ln S_1}{\ln S_{jl}} \pabl{\ln S_{jl}}{S_{jl}}
     +  \pabl{\Delta I_\nu}{S_2} \pabl{S_2}{\ln S_2}
        \pabl{\ln S_2}{\ln S_{jl}} \pabl{\ln S_{jl}}{S_{jl}} \\
    &=& \pabl{\Delta I_\nu}{S_1} S_1 f^{(1)}_{jl} \frac{1}{S_{jl}}
     +  \pabl{\Delta I_\nu}{S_2} S_2 f^{(2)}_{jl} \frac{1}{S_{jl}} \\
    &=& {\rm e}^{-\tau_\nu(s)} \zeta \kappa_1 \Delta s \frac{S_1}{S_{jl}}
        \left(f^{(1)}_{jl}\,g_1(t) + f^{(2)}_{jl}\,g_2(t) \right)
    \label{LambdaCoeff_nonlinear}
\end{eqnarray}
Equation~(\ref{LambdaCoeff_nonlinear}) is used to fill the
$\Lambda$-operator coefficients during a RT formal solution in an
analogous way as described by Eq.~(\ref{LamCoeff}).
Equation~(\ref{LambdaCoeff_nonlinear}) also demonstrates that the
$\Lambda$-operator coefficients $\Lambda_{ik}^{jl}$ now depend on the
source functions.



\begin{figure}
\begin{tabular}{cc}
  \Large $\Lambda$-iteration & \Large $\Lambda$-operator\\
  \hspace*{-11mm}\includegraphics[width=85mm]{epsfig/slab_OLD_conv.eps} &
  \hspace*{-11mm}\includegraphics[width=85mm]{epsfig/slab_LOP_conv.eps} \\
  \hspace*{-11mm}\includegraphics[width=85mm]{epsfig/slab_OLD_resu.eps} &
  \hspace*{-11mm}\includegraphics[width=85mm]{epsfig/slab_LOP_resu.eps} \\
\end{tabular}
\caption{Semi-infinite slab with isotropic irradiation $I_\nu^0$ and
  negligible dust emission ($T_{\rm dust}\!=\!20\,$K). The red dashed
  line is $\frac{1}{2}\exp(-\tau)$, the black line is
  $\frac{1}{2}E_2(\tau)$, where $E_2$ is the second exponential
  integral, the analytic solution of the problem for albedo $a=0$.
  Due to scattering, much more radiation is transported into the
  slab. The cyan line is an approximate solution of the problem,
  $\frac{1}{2} \alpha_1 E_2(\alpha_2\tau)$ with $\alpha_1=1+0.7\cdot
  a$ and $\alpha_2=1-0.25\cdot a -0.65\cdot a^4$, used as initial
  guess. The solution is plotted as thick grey line.}
\label{fig2}
\end{figure}

\begin{figure}
\hspace*{-8mm}\resizebox{18cm}{!}{
\begin{tabular}{ccc}
  $\Lambda$-iteration (proper) 
 &$\Lambda$-iteration (tweaked) 
 &$\Lambda$-operator\\
  \hspace*{-11mm}\includegraphics[width=68mm]{epsfig/RT2_conv_OLDproper.eps} 
 &\hspace*{ -7mm}\includegraphics[width=68mm]{epsfig/RT2_conv_OLDtweak.eps} 
 &\hspace*{ -7mm}\includegraphics[width=68mm]{epsfig/RT2_conv_LOP.eps} 
 \\[-2mm]
  \hspace*{-11mm}\includegraphics[width=68mm,clip]{epsfig/RT2_Jnu_OLDproper.eps}
 &\hspace*{ -7mm}\includegraphics[width=68mm,clip]{epsfig/RT2_Jnu_OLDtweak.eps}
 &\hspace*{ -7mm}\includegraphics[width=68mm,clip]{epsfig/RT2_Jnu_LOP.eps}
 \\[-2mm]
  \hspace*{-11mm}\includegraphics[width=68mm]{epsfig/RT2_Td_OLDproper.eps}
 &\hspace*{ -7mm}\includegraphics[width=68mm]{epsfig/RT2_Td_OLDtweak.eps}
 &\hspace*{ -7mm}\includegraphics[width=68mm]{epsfig/RT2_Td_LOP.eps}
\end{tabular}}
\hspace*{-6mm}\parbox{17.5cm}{\caption{Disk radiative transfer benchmark
    (Pinte et al.\ 2009) for $\tau_{1\mu m}\!=\!10^4$. {\bf The
      l.h.s.~model} shows the old $\Lambda$-iteration method without
    central opacity reduction. This model did not converge, there is
    an obvious problem with the midplane temperatures close to the
    star where the disk is extremely optically thick ($\tau_{\rm
      vert}\approx 1000$). The ``resulting'' temperatures in
    this region are still about the same as the initial guesses. {\bf
      The middle model} used the old $\Lambda$-iteration method with
    reduced central opacities, the standard ProDiMo RT. This model did
    converge and gives about the right temperature structure, but the
    attenuation of $\sim 300-700$\,nm light is too little toward the
    midplane, due to the central opacity reductions. {\bf The
      r.h.s.~model} is the new $\Lambda$ operator method, using both
    temperature and scattering corrections. It results in the proper
    temperatures and mean intensities, but the convergence behavior is
    still not satisfying.}}
\label{fig3}
\end{figure}


\end{document}





