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