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

\setlength{\textheight}{24cm}
\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\bruch#1#2{\frac{#1}{#2}}
    \def\abl#1#2{\bruch {d #1}{d #2}}
    \def\pabl#1#2{\bruch {\partial #1}{\partial #2}}
    \def\plus{${\rm \hspace*{0.9mm}\&\hspace*{0.9mm}}$} 
    \def\etal{${\rm \hspace*{1.0mm}et\,al.\,}$}
    \def\nn{\vec{n}}
    \def\rr{\vec{r}}
    \def\r0{\vec{r}_0}
    \def\kabs{\kappa_\nu^{\rm abs}}
    \def\ksca{\kappa_\nu^{\rm sca}}
    \def\kext{\kappa_\nu^{\rm ext}}
    \def\eg{e.\,g.\ }
    \def\ie{i.\,e.\ }
    \def\nH{n_{\rm\langle H \rangle}}
    \def\AA{\stackrel{\hspace{-0.15ex}\scriptstyle_\circ}
                     {\rm\textstyle A}}
    \def\la{\mathrel{\mathchoice
    {\vcenter{\offinterlineskip\halign{\hfil
    $\displaystyle##$\hfil\cr<\cr\sim\cr}}}
    {\vcenter{\offinterlineskip\halign{\hfil$\textstyle##$\hfil\cr
    <\cr\sim\cr}}}
    {\vcenter{\offinterlineskip\halign{\hfil$\scriptstyle##$\hfil\cr
    <\cr\sim\cr}}}
    {\vcenter{\offinterlineskip\halign{\hfil$\scriptscriptstyle##$\hfil\cr
    <\cr\sim\cr}}}}}
    \def\ga{\mathrel{\mathchoice
    {\vcenter{\offinterlineskip\halign{\hfil
    $\displaystyle##$\hfil\cr>\cr\sim\cr}}}
    {\vcenter{\offinterlineskip\halign{\hfil$\textstyle##$\hfil\cr
    >\cr\sim\cr}}}
    {\vcenter{\offinterlineskip\halign{\hfil$\scriptstyle##$\hfil\cr
    >\cr\sim\cr}}}
    {\vcenter{\offinterlineskip\halign{\hfil$\scriptscriptstyle##$\hfil\cr
    >\cr\sim\cr}}}}}

\begin{document}

\sloppy

\centerline{\huge\bf Continuum Radiative Transfer in ProDiMo}
{\ }\\[-1.7ex]
\centerline{\large Peter Woitke, July 2008}

\section{Definition of Rays}

Let $\r0=(x_0, y_0, z_0)$ denote a point of interest where the mean
intensity $J_\nu(\r0)$ is to be calculated.  The direction of a ray
through $\r0$ is specified by a unit vector pointing backward along
that ray 
\begin{equation}
    \left(\begin{array}{c} n_1 \\ n_2 \\ n_3 \end{array}\right)
  = \left(\begin{array}{c} \sin\theta\cos\phi \\ 
                     \sin\theta\sin\phi \\
                     \cos\theta              \end{array}\right)
\end{equation}
In ProDiMo, there is one ray responsible for the solid angle occupied
by the star and $I\times J$ rays for the remainder of the $4\pi$ solid angle
defined by a 2D-mesh of angular grid points $\{\theta_i\,|\,i=1,...,I\}$
and $\{\phi_j\,|\,j=1,...,J\}$
\begin{eqnarray}
  \theta_i &=& \theta_\star + (\pi-\theta_\star)
               \left(\frac{i-1}{I-1}\right)^{p} \\
  \phi_j   &=& 2\pi \frac{j-1}{J-1} \ ,
\end{eqnarray}
where $\theta_\star=\arcsin(R_\star/r)$ is the half angular diameter
of the star as seen from point $\r0$, $R_\star$ the stellar radius,
and $r=\sqrt{x_0^2+y_0^2+z_0^2}$ the radial distance to the star. The
power index $p>1$ ($p\approx1.5$) assures that there are more rays toward
the hot inner regions than toward the cooler interstellar side. The
integration over solid angle is carried out as
\begin{equation}
  4\pi = \Omega_\star + \sum_{i=1}^{I-1}\sum_{j=1}^{J-1} d\Omega_{ij} \ ,
\end{equation}
where
\begin{eqnarray}
  \Omega_\star &=& 2\pi(1-\cos\theta_\star)\\
  d\Omega_{ij}  &=& (\phi_{j+1}-\phi_j)(\cos\theta_i-\cos\theta_{i+1})
\end{eqnarray}
The direction toward the centre of a solid angle interval
$d\Omega_{ij}$ is given by $\bar{\theta_i}=(\theta_{i+1}+\theta_i)/2$
and $\bar{\phi_j}=(\phi_{j+1}+\phi_j)/2$. In order to achieve that the
direction $(n_1=0, n_2=0, n_3=1)$, \ie $\theta=0$, is pointing toward
the centre of the star, we apply the following turning matrix
\begin{equation}
  \left(\begin{array}{c} n_x \\ n_y \\ n_z \end{array}\right)
   = \left(\begin{array}{ccc}  
            \cos\alpha & 0 & -\sin\alpha \\
                 0     & 1 &      0      \\
            \sin\alpha & 0 & \cos\alpha  \end{array}\right)
     \left(\begin{array}{c} n_1 \\ n_2 \\ n_3 \end{array}\right)
\end{equation}
where $\alpha=\beta+90^o$ and $\tan\beta=z_0/\sqrt{x_0^2+y_0^2}$.

\section{Solution of the Radiative Transfer Equation}

Along each ray in direction $\nn=(n_x,n_y,n_z)$ (backward to the
photon propagation direction), we solve the radiative transfer
equation
\begin{equation}
  \abl{I_\nu}{\tau_\nu} = S_\nu - I_\nu
\end{equation}
assuming LTE and coherent isotropic scattering
\begin{equation}
  S_\nu = \frac{\kabs B_\nu(T)+ \ksca J_\nu}{\kext} \ .
\end{equation}
$I_\nu$ is the spectral intensity, $J_\nu=\frac{1}{4\pi}\int I_\nu
d\Omega$ the mean intensity, $S_\nu$ the source function, and $\kabs$,
$\ksca$ and $\kext=\kabs+\ksca$ are the (dust) absorption, scattering
and extinction coefficients $[\rm cm^{-1}]$. The optical depth 
backward along the ray is given by
\begin{equation}
  \tau_\nu = \int_0^s \kext(\r0+s'\nn)\,ds'
\end{equation}
The formal solution of the transfer equation is
\begin{equation}
  I_\nu = I^{\rm inc}_\nu \exp(-\tau_\nu) 
        + \int_0^{\tau_\nu} S_\nu(\tau_\nu')\exp(-\tau_\nu')\,d\tau_\nu'
  \label{formsol}
\end{equation} 
where $I^{\rm inc}$ is the incident intensity at the end of the
(backward) ray at optical depth $\tau_\nu$. 

Numerically, the integration of the transfer equation is carried out
as follows. We start the ray at $s=0$ with $\tau_\nu=0$ and $I_\nu=0$,
and subdivide the ray into suitable spatial steps $\Delta s$ of the
order of the spatial grid resolution. For each step from
$\rr=\r0+s\,\nn$ to $\r0+(s+\Delta s)\nn$, the optical depths are
incremented for all frequencies $\nu$ according to Simpson's rule
\begin{equation}
  \Delta\tau_\nu = \left(
     \kext\big(\rr\big)
   +4\kext\big(\rr+\frac{1}{2}\Delta s\,\nn\big)
    +\kext\big(\rr+\Delta s\,\nn\big)\right)
  \frac{\Delta s}{6}
\end{equation}
$I_\nu$ is incremented according to Eq.~(\ref{formsol}) as
\begin{equation}
  I_\nu(\tau_\nu+\Delta\tau_\nu) = I_\nu(\tau_\nu)
  + \exp(-\tau_\nu) \int_0^{\Delta\tau_\nu} S_\nu(\tau_\nu+t)\exp(-t)\,dt
\end{equation}
where we assume $S_\nu(\tau_\nu+t)\approx a+bt$ with
$a=S_\nu(\rr)$ and $b=(S_\nu(\rr+\Delta s\,\nn)-a)/\Delta\tau_\nu$, resulting
in
\begin{equation}
  \int_0^{\Delta\tau_\nu} S_\nu(\tau_\nu+t)\exp(-t)\,dt
  \;=\; (a+b) - \rm{e}^{-\Delta\tau_\nu}(a+b+b\Delta\tau_\nu) \ .
\end{equation}
At the end of the ray, where it leaves the model volume, we have to
add the attenuated incident intensity $I^{\rm inc}_\nu
\exp(-\tau_\nu)$ according to Eq.~(\ref{formsol}), where for the core
ray $I^{\rm inc}_\nu=I_\nu^\star$, and for all other rays $I^{\rm
inc}_\nu=I_\nu^{\rm ISM}$. Note that non-core rays may temporarily
leave the disk, but re-enter the disk after some large distance. These
cases -- called ``passages'' in the code -- are treated with large
$\Delta s$ and zero opacity.

The functions $\kext(\rr)$ and $S_\nu(\rr)$ are log-interpolated
between the spatial grid points of ProDiMo. Since axisymmetry is
assumed, a 2D-interpolation is performed in cylinder coordinates
$(\sqrt{x^2+y^2},|z|)$.

\section{The scattering problem}

Since the mean intensities $J_\nu$ enter into the source-functions
which, via the transfer equation, determine the intensities and hence
the mean intensities, a closed solution of the radiative transfer
problem with scattering is not possible, even if the (dust)
temperature structure $T(\rr)$ is known. 

Thus, an iteration is required. During one iteration step in ProDiMo,
the mean intensities are fixed, and the source functions are pre-calculated.
\begin{equation}
  S_\nu(\rr) = \frac{\kabs(\rr) B_\nu\big(T(\rr)\big)
                   + \ksca(\rr) J_\nu^{\rm old}(\rr)}{\kext(\rr)} \ .
  \label{sourcefunc}
\end{equation}
After having solved all rays from all points for all frequencies, the mean
intensities are renewed as
\begin{equation}
  J_\nu(\r0) = \frac{1}{4\pi}\left( I_\nu(\r0,0,0) \Omega_\star + 
               \sum_{i=1}^{I-1}\sum_{j=1}^{J-1} I_\nu(\r0,\theta_i,\phi_j)
               \,d\Omega_{ij} \right)\ ,
\end{equation}
If the maximum change $|\log_{10} J_\nu(\r0)-\log_{10} J_\nu^{\rm
old}(\r0)|$ (all points, all frequencies) is still larger than some
threshold (\eg 0.02), the source functions are updated according to
Eq.~(\ref{sourcefunc}) and the radiative transfer is calculated 
again.

\section{Spectral bands and band-mean quantities}

The purpose of performing a continuum radiative transfer in ProDiMo,
in the end, will be to calculate certain integrals over frequency space,
\eg solving the condition of radiative equilibrium for the dust grains
to obtain $T(r)$, calculating the strength of the UV radiation field
$\chi$ between $912\AA$ and $2050\AA$ etc. The evaluation of these
frequency integrals, in principle, requires a large number of
frequency grid points $\nu_k$ which is computationally expensive.

However, the dust opacities are generally quite smooth and the
solution of the radiative transfer equation on a large frequency grid
will simply yield very similar results for neighbouring points, which
means redundant work. Therefore, it is advantageous to ``interchange''
the order of integration and radiative transfer, and to switch from
a monochromatic treatment to a treatment with spectral bands.

Consider a coarse grid of frequency points $\{\nu_k\,|\,k=0,...,K\}$
(\eg K=10) which covers the whole SED, typically from $1000\AA$ to
$500\mu$m. Instead of $B_\nu(T)$ we consider band means as
\begin{equation}
  B_k(T) = \int_{\nu_{k-1}}^{\nu_k} B_\nu(T)\,d\nu / \Delta\nu_k
\end{equation}
where $\Delta\nu_k=\nu_k-\nu_{k-1}$. In a similar way, we treat the
intensities, mean intensities and opacities
\begin{eqnarray}
  I_k &=& \int_{\nu_{k-1}}^{\nu_k} I_\nu\,d\nu / \Delta\nu_k \\
  J_k &=& \int_{\nu_{k-1}}^{\nu_k} J_\nu\,d\nu / \Delta\nu_k \\
  \kappa_\nu &=& \int_{\nu_{k-1}}^{\nu_k} \kappa_\nu\,d\nu / \Delta\nu_k \\
\end{eqnarray}
Now, we exchange the index $\nu$ in all equations presented so far by
the index $k$, and retrieve the recipes for the band-averaged continuum
radiative transfer. This is of course not an exact treatment, because it
ignores all non-linear couplings, but an approximation that allows us to
use fewer frequency grid points without loosing too much accuracy.
For example, the approximation for the total thermal dust emission 
\begin{equation}
  4\pi \int_{\nu_0}^{\nu_K} \kabs B_\nu(T)\;d\nu 
  \;\approx\; 4\pi \sum_{k=1}^{K} \kappa_k^{\rm abs} B_k(T)\,\Delta\nu_k \\
\end{equation}
is actually exact if $\kabs$ is constant, and still provides a good
approximation if $\kabs$ varies only little within each spectral band.


\section{Irradiation}

If the condition of radiative equilibrium for the dust grains is
accounted for, the radiation field in (and around) the disk is
completely determined by the stellar and interstellar irradiation and
the geometry of the (dust) opacity structure. Therefore, setting the
irradiation as realistic as possible if of ample importance.

\begin{figure}[t]
\centering
\includegraphics[height=6.5cm]{SunSpectra.eps}
\caption{\footnotesize Input stellar spectrum for ProDiMo (red line).
         The example shows a solar model spectrum ($T_{\rm
         eff}=5800\,$K, $\log(g)=4.5$, $Z=1$, black line) and with
         added chromospheric flux from HD 143006 (citation 20xx). Note that
         the ionising and photodissociating radiation is assumed to be
         bracketed by $91.2\,$nm and $205\,$nm.}
\end{figure}

\subsection{Stellar Irradiation}


ProDiMo treats the irradiation from the central star simply by one
core ray for each considered grid point in the model volume which is
assumed to be representative for the solid angle occupied by the
stellar disk (no limb-darkening).  The incident intensity originating
from the stellar atmosphere is deduced from a (modified) model
spectrum from stellar atmosphere codes, \eg a Kurucz-model or a
PHOENIX-model. According to the neglection of limb-darkening, the
incident stellar intensities are related to the surface flux at the
stellar radius via
\begin{equation}
  I_\nu^\star = \bruch{1}{\pi} F_\nu^\star(R_\star) 
              = 4 H_\nu^\star(R_\star) 
\end{equation}
where $F_\nu$ is the flux and $H_\nu$ the Eddington flux. Early-type
stars, which are active and possibly accreting, usually have excess UV
flux as compared to model atmospheres, in particular for cool
stars. This is of central importance for ProDiMo, because it is just
this radiation that ionises and photo-dissociates the atoms and
molecules in the disk. Therefore, we add extra UV-flux from
appropriate sources (observations, recipies ...) to the stellar input
spectrum, see Fig.~1.

The stellar atmosphere fluxes as well as the extra UV can be strongly
peaked and varying in frequency space, especially in the blue. As this
radiation interacts with almost completely flat dust opacities in the
disk, it really makes sense in integrate first (to collect the power in
the spectral bands, see Sect.~4) and then do the radiative transfer
with band mean opacities\footnote{This will be more complicated, 
if gas opacities are important for the continuum radiative transfer, 
see Woitke (2008, MC-paper).}
\begin{equation}
  \mbox{core-rays:}\qquad I_k^{\rm inc} = 
    \int_{\nu_{k-1}}^{\nu_k} I_\nu^\star\,d\nu / \Delta\nu_k
\end{equation}


\subsection{Interstellar Irradiation}

Compared to the stellar irradiation, the interstellar irradiation is
much less meaningfull for the model. Assuming an isotropic
interstellar radiation field, all incident intensities for
non-core rays are approximated by a highly diluted 20000\,K-black-body
field plus the 2.7\,K-cosmic background.
\begin{equation}
  \mbox{non-core-rays:}\qquad I_k^{\rm inc} 
        = B_k(2.7\,{\rm K}) \;+\; \chi^{\rm ISM}\!\cdot 1.71\cdot 
         \mbox{9.85357\,E-17}\cdot B_k(20000\,{\rm K}) 
\end{equation}
The applied dilution factor is close to the value given by
(Draine\plus Bertoldi 1996: 1.06\,E-16). The setting of the dilution
factor for ProDiMo and the additional factor 1.71 will become clearer
in paragraph 6.1.


\section{Applications}

\subsection{$\chi$ from UV radiative transfer}

$\chi$ is a measure of the strength of the UV radiation field, which
directly enters into the calculation of photo-rates in the chemistry.
Draine\plus Bertoldi (1996) define it as
\begin{equation}
  \widetilde{\chi} = \frac{\lambda u_\lambda(1000\AA)}
                          {4\cdot 10^{-14} {\rm erg\,cm^{-3}}}
  \label{chidef1}
\end{equation}
$\lambda$ is the wavelength and $u_\lambda=\frac{4\pi}{c}J_\lambda$
the photon energy density.  This is quite an unfortunate historical
definition, because the photo-cross-sections of atoms and molecules
are by no means restricted to a narrow interval around $1000\,\AA$,
and the UV flux can be strongly varying with wavelength. In
particular, the Ly\,$\alpha$ line at $\lambda=1216\AA$ can contain
most of the UV photon energy. It is hence preferable to work with an
integral definition, and we have to adopt the definition used in the
UMIST chemical reaction database, because we will calculate these
rates with the $\chi$ determined from our radiative
transfer. Unfortunately, Woodall\etal (2007) give no further details about
this issue in their UMIST\,2006 paper. Draine\plus Bertoldi (1996)
write
\begin{equation}
  \widetilde{F} = \frac{1}{h}\int_{912\AA}^{1110\AA} 
                  \lambda u_\lambda\,d\lambda
\end{equation}
which yields a value of $\widetilde{F}=1.22199\,$E+7$\rm\,cm^{-2}\,s^{-1}$ for
Habing's UV field:
\begin{equation}
  \lambda u_\lambda^{\rm Habing} = 10^{-14}{\rm cm^{-2}\,s^{-1}}
  \left(-\frac{25}{6}\lambda_3^3 
        +\frac{25}{2}\lambda_3^2 
        -\frac{13}{3}\lambda_3\right)
\end{equation}
with $\lambda_3=\lambda/1000\AA$. R{\"o}llig\etal (2007) relate
$\chi=1$ to a ``unit Draine field'' and we will follow this idea in
ProDiMo. From the orininal work by Draine (1978), Draine\plus Bertoldi
(1996) deduced
\begin{equation}
   \lambda u_\lambda^{\rm Draine} = 1.71 \times 10^{-14}{\rm cm^{-2}\,s^{-1}}
         \;\frac{31.016\lambda_3^2 - 49.913\lambda_3 + 19.897}{\lambda_3^5}
\end{equation}
The ``magic factor'' 1.71 is because $\lambda u_\lambda^{\rm
Draine}(1000\AA) \approx 1.71 \times \lambda u_\lambda^{\rm
Habing}(1000\AA) = 1.71 \times 10^{-14}{\rm cm^{-2}\,s^{-1}}$, as can
be seen from Fig.~3 of Draine (1978). Consequently, $\chi$ should be
equal to 1.71 according to the historic definition (Eq.~\ref{chidef1})
and ``a unit Draine field'' would imply to use $\chi=1.71$ for the
UMIST photo-rates as suggested by Draine\plus Bertoldi
(1996). However, we doubt that this is correct and will put $\chi=1$
for a unit Draine field. Acordingly, our (working) definition for
$\chi$ is
\begin{equation}
  \chi = F/F_{\rm Draine} 
\end{equation}
with
\begin{eqnarray}
  F &=& \frac{1}{h}\int_{912\AA}^{2050\AA} \lambda u_\lambda\,d\lambda \\
  F_{\rm Draine} &=& 1.921457\,\mbox{E+8}\rm\,cm^{-2}\,s^{-1}\ .
\end{eqnarray}
It is important to note that ``$F$'' is not a flux (a flux would be
zero in an isotropic medium!) but a measure related to the mean
intensity of the radiation field. Given the definition for
$u_\nu=\frac{4\pi}{c}J_\nu$ is it straightforward to show that
\begin{equation}
  \chi = \frac{4\pi}{h\bar{\nu}_1} \,J_1 \Delta\nu_1 \,\Big{/}\, F_{\rm Draine}
\end{equation}
with $\bar{\nu}_k=\sqrt{\nu_k\nu_{k-1}}$, because we choose band
1 to be bracketed by $912\AA$ and $2050\AA$.



\subsection{Dust temperature from radiative transfer}

The condition of radiative equilibrium for the dust component, with the
applied approximation of band means, is given by
\begin{equation}
  \sum \kappa_k^{\rm abs} J_k(\rr)    \Delta\nu_k =
  \sum \kappa_k^{\rm abs} B_k\big(T(\rr)\big) \Delta\nu_k \ ,
\end{equation}
which states an implicit equation for the determination of the dust
temperature $T(\rr)$, which is solved by a Newton-Raphson iteration.

The inclusion of the determination of the dust temperature structure
is optional in ProDiMo and slows down the convergence of the radiative
transfer considerably.
\begin{equation}
  J_k(\rr) \begin{array}{c} \to T(\rr) \to\\
           {\resizebox{1.9cm}{!}{\mbox{$\longrightarrow$}}}
           \end{array} 
  S_k(\rr) \stackrel{\rm RT}{\longrightarrow} J_k(\rr) 
\end{equation} 
 In order to accelerate the convergence, we apply the procedure of
 Auer (1984) for this $\lambda$-type iteration scheme.



\subsection{Backround mean intensities for non-LTE calculations}

The non-LTE modelling of atoms, ions and molecules depends not only on
the kinetic gas temperature, the particle densities and the collision
partner densities, but also on the continuum mean intensities at the
wavelengths of the spectral lines, which originates in general from
distant sources like the central star, the ISM, or the warm dust in
the disk.  The ``background mean intensities'' $J_\nu^{\rm back}$ are
identified to be just given by the local mean intensities calculated
from the dust continuum radiative transfer persented here. For strong
radiations fields and cool temperatures, radiative line cooling even
turns into radiative heating due to line absorption followed by
collisional de-excitation. Thus, setting the background intensities
properly has an important impact on the gas energy balance.

\begin{figure}[htbp]
\centering
\includegraphics[height=7.5cm,width=9cm]{Jcalc.eps}
\caption{\footnotesize Calculated band-mean mean intensities at one
         point in a model (12 black dots) and a cubic spline interpolation
         through these points (black line). The vertical lines
         indicate the interval boundaries of the 12 spectral bands
         considered. The red line shows the band-mean incident stellar
         intensities $\nu I_\nu^\star\,\rm[erg\,cm^{-2}\,s^{-1}\,sr^{-1}]$
         and the blue line shows the incident interstellar intensities. The
         radiation field has basically two components, the dust
         attenuated UV -- near IR part, orgininating mainly from the
         star, and the self-generated mid -- far IR part, origiating 
         from the thermal dust emission.}
\end{figure}

Thus, we can ``feed'' our non-LTE models with the calculated
$J_\nu(\rr)$. In order to obtain the required monochromatic mean
intensities at the line center positions without loosing too much
accurancy, we apply a local cubic spline interpolation in frequency
space as depicted in Fig.~2.


\end{document}
