\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}