Free Access
Volume 501, Number 1, July I 2009
Page(s) 383 - 406
Section Numerical methods and codes
Published online 29 April 2009

Radiation thermo-chemical models of protoplanetary disks

I. Hydrostatic disk structure and inner rim

P. Woitke1,2 - I. Kamp3 - W.-F. Thi4

1 - UK Astronomy Technology Centre, Royal Observatory, Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK
2 - School of Physics & Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK
3 - Kapteyn Astronomical Institute, Postbus 800, 9700 AV Groningen, The Netherlands
4 - SUPA[*], Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

Received 10 February 2009 / Accepted 1 April 2009

Context. Emission lines from protoplanetary disks originate mainly in the irradiated surface layers, where the gas is generally warmer than the dust. Therefore, interpreting emission lines requires detailed thermo-chemical models, which are essential to converting line observations into understanding disk physics.
Aims. We aim at hydrostatic disk models that are valid from 0.1 AU to 1000 AU to interpret gas emission lines from UV to sub-mm. In particular, our interest lies in interpreting far IR gas emission lines, such as will be observed by the Herschel observatory, related to the G ASPS open time key program. This paper introduces a new disk code called P ROD IM O.
Methods. We combine frequency-dependent 2D dust continuum radiative transfer, kinetic gas-phase and UV photo-chemistry, ice formation, and detailed non-LTE heating & cooling with the consistent calculation of the hydrostatic disk structure. We include Fe  II and CO ro-vibrational line heating/cooling relevant to the high-density gas close to the star, and apply a modified escape-probability treatment. The models are characterised by a high degree of consistency between the various physical, chemical, and radiative processes, where the mutual feedbacks are solved iteratively.
Results. In application to a T Tauri disk extending from 0.5 AU to 500 AU, the models show that the dense, shielded and cold midplane ( $z/r \la 0.1$, $T_{\hspace*{-0.2ex}\rm g}\approx T_{\hspace*{-0.2ex}\rm d}$) is surrounded by a layer of hot ( $T_{\hspace*{-0.2ex}\rm g}\approx 5000~$K) and thin ( $n_{{\rm\langle H\rangle}}\!\approx\!10^{~7}$ to $10^{~8}\rm ~cm^{-3}$) atomic gas that extends radially to about 10 AU and vertically up to $z/r\!\approx\!0.5$. This layer is predominantly heated by the stellar UV (e.g. PAH-heating) and cools via Fe  II semi-forbidden and O I 630 nm optical line emission. The dust grains in this ``halo'' scatter the starlight back onto the disk, which affects the photochemistry. The more distant regions are characterised by a cooler flaring structure. Beyond $r \ga 100~$AU, $T_{\hspace*{-0.2ex}\rm g}$ decouples from  $T_{\hspace*{-0.2ex}\rm d}$ even in the midplane and reaches values of about $T_{\hspace*{-0.2ex}\rm g}\approx 2T_{\hspace*{-0.2ex}\rm d}$.
Conclusions. Our models show that the gas energy balance is the key to understanding the vertical disk structure. Models calculated with the assumption $T_{\hspace*{-0.2ex}\rm g}\!=\!T_{\hspace*{-0.2ex}\rm d}$ show a much flatter disk structure. The conditions in the close regions (<10 AU) with densities $n_{{\rm\langle H\rangle}}\!\approx\!10^{~8}$ to $10^{~15}\rm ~cm^{-3}$ resemble those of cool stellar atmospheres and, thus, the heating and cooling is more like in stellar atmospheres. The application of heating and cooling rates known from PDR and interstellar cloud research alone can be misleading here, so more work needs to be invested to identify the leading heating and cooling processes.

Key words: astrochemistry - radiative transfer - methods: numerical - stars: formation - stars: circumstellar matter

1 Introduction

The structure and composition of protoplanetary disks play a key role in understanding the process of planet formation. From thermal and scattered light observations, we know that protoplanetary disks are ubiquitous in star-forming regions and that the dust in these disks evolves on timescales of 106 yr (Watson et al. 2007; Haisch et al. 2006). However, dust grains represent only about 1% of the mass in these disks - 99% of their mass is gas.

While observations of the dust in these systems have a long history (e.g. Beckwith et al. 1990), gas observations are intrinsically more difficult and until recently have focused on rotational lines of abundant molecules such as CO, HCN, HCO+(e.g. Thi et al. 2004; Beckwith et al. 1986; Dutrey et al. 1997; van Zadelhoff et al. 2001; Koerner et al. 1993). These lines generally trace the outer cooler regions of protoplanetary disks and probe a layer at intermediate heights, where the stellar UV radiation is sufficiently shielded to suppress photodissociation, but still provides enough ionisation to drive a rich ion-molecule chemistry (Bergin et al. 2007). However, the sensitivity of current radio telescopes only allows observations of small samples (Dent et al. 2005) and, in a few cases, detailed studies of individual objects (e.g. Qi et al. 2003; Semenov et al. 2005; Qi et al. 2008). The gas temperature in the disk surface down to continuum optical depth of $\sim$1decouples from the dust temperature and ranges from a few thousand to a few hundred Kelvin (Jonkheid et al. 2004; Kamp & Dullemond 2004; Dullemond et al. 2007). The hot inner disk can show ro-vibrational line emission in the near-IR either due to fluorescence (at the disk surface) and/or thermal excitation (e.g. Bitner et al. 2007; Bary et al. 2003; Brittain et al. 2007). More recently, near-IR gas lines have also been detected in Spitzer IRS spectra, revealing the presence of water, H2 and the importance of X-rays (Lahuis et al. 2007; Salyk et al. 2008; Pascucci et al. 2007). The launch of the Herschel satellite in 2009, opens yet another window to study the gas component of protoplanetary disks through the dominant cooling lines [O  I], [C II] at the disk surface, as well as many additional molecular tracers of the warmer inner disk such as water and CO. The study of the gas in protoplanetary disks is the main topic of the Herschel open time Key Program ``Gas in Protoplanetary Systems'' (G ASPS, PI: Dent). Other guaranteed and open time Key Programs, such as ``Water in Star Forming Regions with Herschel'' (W ISH, PI: van Dishoeck) and ``HIFI Spectral Surveys of Star Forming Regions'' (PI: Ceccarelli), will also observe gas lines in a few disks.

Disk structure modelling was initially driven by dust observations and developed from a simple two-layer disk model (e.g. Chiang & Goldreich 1997) into detailed dust continuum radiative transfer models that are coupled with hydrostatic equilibrium (e.g. Pinte et al. 2006; Dullemond et al. 2002; D'Alessio et al. 1998; Dullemond & Dominik 2004). The assumption in all these models is that gas and dust are well coupled and the hydrostatic scale height then follows from the dust temperature. However, the gas temperature decouples from the dust temperature and the vertical disk structure will adjust to the gas scale height, forcing the dust to follow if it is dynamically coupled. This approach has been followed by Nomura & Millar (2005) and Gorti & Hollenbach (2008,2004). Nomura & Millar use a small chemical network (only CO, C+ and O) and a limited number of heating/cooling processes namely photoelectric heating, gas-grain collisions and line cooling from [O  I], [C II] and CO. Gorti & Hollenbach use an extended set of reactions (84 species, $\sim$600 reactions) and the relevant low-density heating/cooling processes drawn from photo dissociation region (PDR) physics. Other models do not solve for the vertical hydrostatic disk structure (e.g.  Woods & Willacy 2008; Kamp & Dullemond 2004; Meijerink et al. 2008). Kamp & Dullemond use a chemical reaction network of $\sim$250 reactions among 48 species and a set of heating/cooling processes comparable to Gorti & Hollenbach (2004). The models of Meijerink ${\rm\hspace*{0.8ex}et\hspace*{0.7ex}al.\hspace*{0.7ex}}$focus entirely on the X-ray irradiation of the disk, thus excluding UV processes; the chemical reaction network is limited to 25 species and 125 reactions. Woods ${\rm\hspace*{0.7ex}\&\hspace*{0.7ex}}$Willacy use again a standard set of PDR heating/cooling processes, but also account for X-rays. Their chemical network includes 475 gas and ice species connected through $\sim$8000 gas phase and surface reactions.

This paper presents a new disk code that includes additional heating/cooling processes relevant for the high densities and high temperatures present in the inner parts of the disk, resembling the conditions in tenuous atmospheres of cool stars. The models are characterised by a high degree of consistency between the various physical, chemical and radiative processes. In particular, the results of a full 2D dust continuum radiative transfer are used as input for the UV photo-processes and as radiation background for the non-LTE modelling of atoms and molecules to calculate the line heating and cooling rates. This allows the models to extend closer to the star and include modelling of the so-called inner rim.

The paper is structured as follows. Section 2 introduces the new code P ROD IM O and presents the concept of global iterations. Section 3 describes the assumptions used to calculate the hydrostatic disk structure including ``soft edges''. In Sect. 4, we present the 2D dust continuum radiative transfer with scattering and band-mean opacities. Section 5 summarises the gas-phase and photo-chemistry dependent on the UV continuum transfer results. In Sect. 6, we outline the heating and cooling rates included in our model and present a modified escape-probability method. Section 7 closes the theory part of the paper with the calculation of the sound speeds as preparation of the next calculation the the disk structure. We apply P ROD IM O to a standard T Tauri-type protoplanetary disk with disk mass $0.01~M_\odot$ which extends from 0.5 AU to 500 AU in Sect. 8. The resulting physical and chemical structure of the disk is shown and compared to a model where we assume $T_{\hspace*{-0.2ex}\rm g}\!=\!T_{\hspace*{-0.2ex}\rm d}$. We conclude the paper in Sect. 9 with an outlook to future applications.

2 ProDiMo

P ROD IM O is an acronym for Protoplanetary Disk Model. It is based on the thermo-chemical models of Inga Kamp (Kamp & Bertoldi 2000; Kamp & Dullemond 2004; Kamp & van Zadelhoff 2001), but completely re-written to be more flexible and to include more physical processes.

P ROD IM O uses global iterations to consistently calculate the physical, thermal and chemical structure of protoplanetary disks. The iterations involve 2D dust continuum radiative transfer, gas-phase and photo-chemistry, thermal energy balance of the gas, and the calculation of the hydrostatic disk structure in axial symmetry (see Fig. 1). The different components will be explained separately in the forthcoming sections.

Physical processes not yet included are X-ray heating, X-ray chemistry, spatially dependent dust properties, and PAH-chemistry. These processes will be addressed in future papers. P ROD IM O is under current development. The code can be downloaded from, start at to get a P ROD IM O user account.

\end{figure} Figure 1:

Concept of global iterations in P ROD IM O. The circular arrows on the r.h.s. indicate sub-iterations. For example, the dust temperature structure needs to be iterated in the continuum radiative transfer.

Open with DEXTER

3 Hydrostatic disk structure

We consider the hydrostatic equation of motion in axial symmetry with rotation around the z-axis, but $v_r\!=\!0$ and $v_z\!=\!0$

$\displaystyle \frac{v_\phi^2}{r} = \frac{1}{\rho} \frac{\partial p}{\partial r} + \frac{\partial \Phi}{\partial r}$     (1)
$\displaystyle 0 = \frac{1}{\rho} \frac{\partial p}{\partial z} + \frac{\partial \Phi}{\partial z}
,$     (2)

where vr, vz and $v_\phi$ are the three components of the velocity field, p the gas pressure and $\rho$ the mass density, respectively. $r\!=\!(x^2+y^2)^{1/2}$ is the distance from the symmetry axis and z is the distance from the midplane. Neglecting self-gravity, the gravitational potential is given by

\begin{displaymath}\Phi(r,z) = -~\frac{G M_\star}{\sqrt{r^2+z^2}}
\end{displaymath} (3)

where $M_\star$ is the stellar mass and G the gravitational constant. We follow the idea of ``1+1D'' modelling (Malbet et al. 2001; Dullemond et al. 2002; D'Alessio et al. 1998) by assuming that the radial pressure gradient $1/\rho(\partial p/\partial r)$ in Eq. (1) is small compared to centrifugal acceleration and gravity, in which case the radial and the vertical components of the equation of motion decouple from each other. The radial component then simply results in circular Keplerian orbits

\begin{displaymath}v_\phi = \left( \frac{ r^2~GM_\star}{\big(r^2+z^2\big)^{3/2}}
\right)^{1/2} ,
\end{displaymath} (4)

leaving the radial distribution of matter undetermined, as it is in fact mostly determined by the actual distribution of angular momentum in the disk. Consequently, the vertical component of the equation of motion (Eq. (2)) can be solved independently for every vertical column in the disk

\begin{displaymath}\frac{1}{\rho}\frac{{\rm d}p}{{\rm d}z} = -~\frac{z\;GM_\star}{\big(r^2+z^2\big)^{3/2}}
\end{displaymath} (5)

Equation (5) is integrated from the midplane upwards by substituting the density for the pressure via $p\!=\!c_T^2~\rho$, and assuming that the isothermal sound speed cT is a known function of z. Numerically, we perform this integration by means of an ordinary differential equation (ODE) solver, using a simple pointwise linear interpolation of cT2(z) between calculated grid points cT2(rj,zk). Since for known cT(z) the solution p(z) has a free factor, we put $p(0)\!=\!1$ and scale the results later to achieve any desired column density at distance r

\begin{displaymath}\Sigma(r) = 2\!\!\int\limits_0^{z_{\rm max}(r)}\!\!\rho(r,z)~{\rm d}z ,
\end{displaymath} (6)

where the factor 2 is because of the lower half of the disk, which is assumed to be symmetric. In this paper, we assume a powerlaw distribution of the column density

\begin{displaymath}\Sigma(r) = \Sigma_0\;r^{~-\epsilon}
\end{displaymath} (7)

in the main part of the disk, except for the ``soft edges'' (see Sect. 3.1), and determine $\Sigma_0$ from the specified disk mass  $M_{\rm disk}$

\begin{displaymath}M_{\rm disk} = 2\pi \int\limits_{R_{\rm in}}^{R_{\rm out}} \Sigma(r)~r~{\rm d}r ,
\end{displaymath} (8)

where $R_{\rm in}$ is the inner radius and $R_{\rm out}$ the outer radius of the disk. In summary, supposed that cT2(r,z) is known, the disk structure is determined by the parameters $M_{\rm disk}$, $M_\star$, $R_{\rm in}$, $R_{\rm out}$, and $\epsilon$.

3.1 Soft edges

The application of a radial surface density powerlaw (Eq. (7)) in the disk between $R_{\rm in}$ and $R_{\rm out}$ is, although widely used, obviously quite artificial and even unphysical. Equation (1) demonstrates that an abrupt radial cutoff would produce an infinite force because of the radial pressure gradient $\partial p/\partial r$, which would push gas inward at  $R_{\rm in}$, and outward at  $R_{\rm out}$, respectively, causing a smoothing of the radial density structure at the boundaries.

Let us consider an abrupt cutoff in the beginning and study the motion of the gas as it is pushed inward due to the radial pressure gradient at the inner boundary. Since the specific angular momentum $L_3(R_{\rm in})\!=\!R_{\rm in}~v_\phi(R_{\rm in})\!=\!r~v_\phi(r)$ is conserved during this motion, the gas will spin up as it is pushed inward, until the increased centrifugal force balances the radial pressure gradient (+ gravity). According to Eq. (1) the force equilibrium in this relaxed state is given by

\begin{displaymath}\frac{L_3(R_{\rm in})^2}{r^3} = \frac{1}{\rho}\frac{\partial p}{\partial r} + \frac{\partial \Phi}{\partial r}
\end{displaymath} (9)

which provides an equation for the desired density structure $\rho(r)$. Using $p(r)=c_T^2~\rho(r)$ and assuming $c_T^2\!=\!\rm const$, the result is

\begin{displaymath}\ln\frac{\rho(r_0)}{\rho(R_{\rm in})} ~=~
...L_3(R_{\rm in})^2}{2r^2}
+\Phi(r)\right]^{~r_0}_{~R_{\rm in}}
\end{displaymath} (10)

where r0 is an arbitrary point inside $R_{\rm in}$. Generalising Eq. (10) to column densities (with cT2 measured in the midplane) we write

\begin{displaymath}\Sigma(r_0) ~\approx~ \Sigma(R_{\rm in})~\exp\left(-\frac{1}{...
M_\star}{r}\right]^{~r_0}_{~R_{\rm in}}\right).
\end{displaymath} (11)

A similar expression can be found for the column density outside of the outer boundary. The CO observations of Hughes et al. (2008) show that such treatments can improve model fits. However, we have chosen to apply our approach for soft edges only to the inner boundary in this paper.

To summarise, if angular momentum is transported inside-out in the disk, the density structure may decrease more gradually or even increase further inward (Hartmann et al. 1998). However, it is hard to figure out any circumstances where the column density could decrease more rapidly at the inner rim as compared to Eq. (11).

4 Continuum radiative transfer

The chemistry and the heating & cooling balance of the gas in the disk (see Sects. 5 and 6) depend on the local continuous radiation field $J_\nu(r,z)$ and the local dust temperature $T_{\hspace *{-0.2ex}\rm d}(r,z)$ which is a result thereof. These dependencies include

thermal accommodation between gas and dust, which is usually the dominant heating/cooling process for the gas in the midplane ( $\to\!T_{\hspace*{-0.2ex}\rm d}$);
photo-ionisation and photo-dissociation of molecules, as well as heating by absorption of UV photons, e.g. photo-electric heating ( $\to \!J_{\rm UV}$);
radiative pumping of atoms and molecules by continuum radiation which alters the non-LTE population and cooling rates, sometimes turning cooling into heating ( $\to\!J_\nu$);
surface chemistry on grains, in particular the H2-formation, and ice formation and desorption ( $\to\!T_{\hspace*{-0.2ex}\rm d},J_{\rm UV}$).
Previous chemical models have often treated these couplings by means of simplifying assumptions and approximate formula (e.g. Kamp & Bertoldi 2000; Nomura & Millar 2005; Hollenbach et al. 1991). For a rigorous solution, a full 2D continuum radiative transfer must be carried out, which provides $T_{\hspace *{-0.2ex}\rm d}(r,z)$ and $J_\nu(r,z)$, including the UV part, at every location in the disk.

P ROD IM O solves the 2D dust continuum radiative transfer of irradiated disks by means of a simple, ray-based, long-characteristic, accelerated $\Lambda$-iteration method. From each grid point in the disk, a number of rays (typically about 100) are traced backward along the photon propagation direction, while solving the radiative transfer equation

\begin{displaymath}\frac {{\rm d} I_\nu}{{\rm d} \tau_\nu} = S_\nu - I_\nu
\end{displaymath} (12)

assuming LTE and coherent isotropic scattering

\begin{displaymath}S_\nu = \frac{\kappa_\nu^{\rm abs}B_\nu(T_{\hspace*{-0.2ex}\rm d})+ \kappa_\nu^{\rm sca}J_\nu}{\kappa_\nu^{\rm ext}} \cdot
\end{displaymath} (13)

$I_\nu$ is the spectral intensity, $J_\nu\!=\!\frac{1}{4\pi}\int I_\nu
~\rm d\Omega$ the mean intensity, $S_\nu$ the source function, $B_\nu$the Planck function, and $\kappa_\nu^{\rm abs}$, $\kappa_\nu^{\rm sca}$ and $\kappa_\nu^{\rm ext}=\kappa_\nu^{\rm abs}+\kappa_\nu^{\rm sca}$ $[\rm cm^{-1}]$ are the dust absorption, scattering and extinction coefficients, respectively.

The dust grains of various sizes at a certain location in the disk are assumed to have a unique temperature $T_{\hspace*{-0.2ex}\rm d}$ in modified radiative equilibrium

\begin{displaymath}\Gamma_{\rm dust} ~+ \int \kappa_\nu^{\rm abs}J_\nu~{\rm d}\n...
...ppa_\nu^{\rm abs}B_\nu(T_{\hspace*{-0.2ex}\rm d})~{\rm d}\nu ,
\end{displaymath} (14)

where the additional heating rate $\Gamma_{\rm dust}$ accounts for non-radiative heating (negative for cooling) processes like thermal accommodation with gas particles and frictional heating. An accelerated $\Lambda$ scheme is used to get converged results concerning $J_\nu(r,z)$ and $T_{\hspace *{-0.2ex}\rm d}(r,z)$. The details of this method will be described in the following sections.

4.1 Geometry of rays

Let $\vec{r}_0\!=\!(x_0, y_0, z_0)$ denote a point in the disk where the mean intensities $J_\nu(\vec{r}_0)$ are to be calculated. The direction of a ray starting from $\vec{r}_0$ is specified by a unit vector which points in the reverse direction of the photon propagation

\begin{displaymath}\left(\begin{array}{c} n_1 \\ n_2 \\ n_3 \end{array}\right)
...hi \\
\sin\theta\sin\phi \\
\cos\theta \end{array}\right)
\end{displaymath} (15)

as specified in a local coordinate system where (0,0,1) points toward the star. One ray is reserved for the solid angle occupied by the star $\Omega_\star$ as seen from point $\vec{r}_0$, subsequently called the ``core ray''. All other $I\times J$ rays represent the remainder of the $4\pi$ solid angle by a 2D-mesh of angular grid points $\{\theta_i~\vert~i=1,...,I\}$ and $\{\phi_j~\vert~j=1,...,J\}$
                    $\displaystyle \theta_i$ = $\displaystyle \theta_\star + (\pi-\theta_\star)
\left(\frac{i-1}{I-1}\right)^{p}$ (16)
$\displaystyle \phi_j$ = $\displaystyle \pi~ \frac{j-1}{J-1} ,$ (17)

where $\theta_\star\!= \arcsin(R_\star/d)$ is the half angular diameter of the star as seen from point $\vec{r}_0$, $R_\star$ the stellar radius, and $d\!=\!(x_0^2+y_0^2+z_0^2)^{1/2}$ the radial distance. $\phi$ only ranges from 0 to $\pi$, because the disk problem is symmetric $I_\nu(\vec{r}_0,\theta,+\phi)\!=\!I_\nu(\vec{r}_0,\theta,-\phi)$. A power index $p\!>\!1$ ( $p\!\approx\!1.5$) assures that there are more rays pointing toward the hot inner regions than toward the cooler interstellar side. The integration over solid angle is carried out as
                           $\displaystyle 4\pi$ = $\displaystyle \Omega_\star + \sum_{i=1}^{I-1}\sum_{j=1}^{J-1} {\rm d}\Omega_{ij}$ (18)
$\displaystyle \Omega_\star$ = $\displaystyle 2\pi~(1-\cos\theta_\star)$ (19)
$\displaystyle {\rm d}\Omega_{ij}$ = $\displaystyle 2(\phi_{j+1}-\phi_j)~(\cos\theta_i-\cos\theta_{i+1}) .$ (20)

The central direction of 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$, and these are the angles actually considered in Eq. (15). In order let the core ray with $\theta\!=\!0$ point toward the star, we apply the following rotation matrix

\begin{displaymath}\left(\begin{array}{c} n_x \\ n_y \\ n_z \end{array}\right)
\left(\begin{array}{c} n_1 \\ n_2 \\ n_3 \end{array}\right)
\end{displaymath} (21)

where $\alpha=\beta+90^o$ and $\tan\beta=z_0/(x_0^2+y_0^2)^{1/2}$.

4.2 Solution of the radiative transfer equation

From every grid point $\vec{r}_0$ along each ray in direction $\vec{n}\!=\!(n_x,n_y,n_z)$ we solve the radiative transfer equation (Eq. (12)) backward to the photon propagation direction. The optical depth along the ray is given by

\begin{displaymath}\tau_\nu(s) = \int_0^s \kappa_\nu^{\rm ext}(\vec{r}_0+s'\vec{n})~{\rm d}s'.
\end{displaymath} (22)

The formal solution of the transfer equation Eq. (12) is

\begin{displaymath}I_\nu ~=~ I^{\rm inc}_\nu {\rm e}^{-\tau_\nu(s_{\rm max})}
...ppa_\nu^{\rm ext}(s)~S_\nu(s)~{\rm e}^{-\tau_\nu(s)}\;{\rm d}s
\end{displaymath} (23)

where $I^{\rm inc}_\nu$ is the intensity incident from the end of the ray at $s_{\rm max}$. We start a ray at $s\!=\!0$ with $\tau_\nu\!=\!0$ and $I_\nu\!=\!0$ and choose a suitable spatial step size $\Delta s$. For each step, the opacities and source functions (all wavelengths) at the start point and the end point of the step, $\vec{r}_0\!+\!s~\vec{n}$ and $\vec{r}_0\!+\!(s\!+\!\Delta s)\vec{n}$, are interpolated from the pre-calculated values on the grid points using a 2D-interpolation in cylinder coordinates ((x2+y2)1/2,|z|). For the numerical integration of Eqs. (22) and (23), we assume $\kappa_\nu^{\rm ext}\!=\!A+Bs$and $S_\nu\!=\!C\exp(Ds)$, where the coefficients A,B,C,D are determined by the start and end point values. Simplifying the exponent by putting $\bar{\kappa}^{\rm ext}_\nu\!\approx\!A+B\Delta
s/2$ yields
                           $\displaystyle \tau_\nu(s\!+\!\Delta s)$ = $\displaystyle \tau_\nu(s)
+ \!\int_0^{\Delta s}\hspace*{-1mm}(A+Bs')~{\rm d}s'$ (24)
$\displaystyle I_\nu(s\!+\!\Delta s)$ = $\displaystyle I_\nu(s)
+ C~{\rm e}^{-\tau_\nu(s)}\!\!
\int_0^{\Delta s}\hspace*{-1mm} (A+Bs')\;
{\rm e}^{~(D~-~\bar{\kappa}^{\rm ext}_\nu)~s'}~{\rm d}s' .$ (25)

The numerical integration is carried out with analytic expressions for these integrals. The procedure is repeated for two half steps of size $\Delta s/2$. If the results differ too much, the step size $\Delta s$ is reduced and the step is re-calculated. In case of small differences, the step size is increased for the following step.

At the end of each ray, the attenuated incident intensities $I^{\rm
inc}_\nu {\rm e}^{-\tau_\nu(s_{\rm max})}$ are added according to Eq. (23), where for the core ray the stellar intensity $I^{\rm inc}_\nu=I_\nu^\star$ is used, and for all other rays the interstellar intensity $I^{\rm inc}_\nu=I_\nu^{\rm ISM}$ is applied. Non-core rays may temporarily leave the disk, but re-enter the disk after some large distance. These ``passages'' are treated with large, exactly calculated $\Delta s$ and zero opacities.

For the 2D-interpolation, it turned out to be important to use a log-interpolation for the source function $S_\nu(r,z)$ which can change by orders of magnitude, e.g. across a shadow, within one step. In case of linear interpolation, the numerical radiative transfer shows much more numerical diffusion.

4.3 Irradiation

The radiation field in (and around) passive disks 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 is of prime importance.

\end{figure} Figure 2:

Incident stellar intensity compiled from two sources: a P HOENIX solar model spectrum with $T_{\rm
eff}\!=\!5800~$K, $\log~g\!=\!4.5$, $Z\!=\!1$ (black line) and the chromospheric flux of ``young sun'' HD 129 333 (Dorren & Guinan 1994, red line). Note that the ionising and photo-dissociating flux is assumed to be restricted to the interval [91.2 nm,205 nm] which includes Ly$\alpha $ at 121.6 nm.

Open with DEXTER

Stellar irradiation

For the incident stellar irradiation, a model spectrum from stellar atmosphere codes is used, e.g. a P HOENIX-model[*]. Neglecting limb-darkening, the incident stellar intensities are related to the surface flux at the stellar radius via

\begin{displaymath}I_\nu^\star = \frac{1}{\pi} F_\nu^\star(R_\star)
~=~ 4~H_\nu^\star(R_\star)
\end{displaymath} (26)

where $F_\nu^\star$ is the spectral flux and $H_\nu$ the Eddington flux. Young stars, which are active and possibly accreting, have excess UV as compared to model atmospheres, in particular cool stars. This is of central importance for P ROD IM O , because it is just this radiation that ionises and photo-dissociates the atoms and molecules in the disk. Therefore, we add extra UV-flux as e.g. reduced from observations or given by other recipes, see Fig. 2.

Interstellar irradiation

Assuming an isotropic interstellar radiation field, all incident intensities for non-core rays are approximated by a highly diluted 20 000 K-black-body field plus the 2.7 K-cosmic background.

$\displaystyle I_\nu^{\rm ISM}
= \chi^{\rm ISM}\cdot 1.71\cdot W_{\rm dil}~B_\nu(20~000{\rm ~K})
~+ B_\nu(2.7{\rm ~K}).$     (27)

The applied dilution factor $W_{\rm dil}=9.85357\times10^{-17}$ is calculated from the normalisation $\chi\!=\!1$ according to Eq. (41), which is close to the value given by Draine & Bertoldi (1996). $\chi^{\rm ISM}$ is a free parameter which describes the strength of the UV field with respect to standard interstellar conditions.

4.4 Iteration and dust temperature determination

In order to solve the condition of the dust radiative equilibrium (Eq. (14) and the scattering problem, a simple $\Lambda$-type iteration is applied. The source functions are pre-calculated on the grid points according to Eq. (13), with $J_\nu\!=\!J_\nu^{\rm old}$ and $T_{\hspace*{-0.2ex}\rm d}\!=\!T_{\hspace*{-0.2ex}\rm d}^{\rm old}$, and fixed during one iteration step. After having solved all rays from all points for all frequencies, the mean intensities are updated as

\begin{displaymath}J_\nu(\vec{r}_0) = \frac{1}{4\pi}\left( I_\nu(\vec{r}_0,0,0)~...
~d\Omega_{ij} \right) ,
\end{displaymath} (28)

and the dust temperatures are renewed according to Eq. (14). If the maximum relative change $\vert J_\nu-J_\nu^{\rm old}\vert/(J_\nu+J_{\rm
small})$ (all points, all frequencies, $J_{\rm small}\!=\!10^{-30}\rm
erg~ cm^{-2}~ s^{-1}~ Hz^{-1}~ sr^{-1}$) is larger than some threshold ($\sim$0.01), the source functions are re-calculated and the radiative transfer is solved again. In order to accelerate the convergence, we apply the procedure of Auer (1984). We benchmarked results of our radiative transfer method against results of other Monte-Carlo and ray-based methods in Pinte et al. (2009). The convergence in optically thick disks is tricky, but we can manage test problems up to a midplane optical depths of about $\tau\!=\!10^5$ with this code (see Fig. 3).

\end{figure} Figure 3:

Benchmark for the dust continuum radiative transfer part. Vertical cuts of the calculated dust temperature structure $T_{\hspace *{-0.2ex}\rm d}(r,z)$ are shown at different radii of a disk with midplane optical depth $\tau_{\rm 1\mu{\rm m}}\!=\!10^5$ at $1~\mu$m. The P ROD IM O results are shown with blue diamonds in comparison to the MCFOST results (Pinte et al. 2006) with black lines.

Open with DEXTER

4.5 Spectral bands and band-mean quantities

The main purpose of the continuum radiative transfer in P ROD IM O is to calculate certain frequency integrals, e.g. solving the condition of radiative equilibrium for the dust grains (Eq. (12)) or calculating the local strength of the UV radiation field $\chi $ (Eq. (41)). The incident stellar spectrum is strongly varying in frequency space, especially in the blue and UV (see Fig. 2) and the evaluation of these integrals, in principle, requires a large number of frequency grid points $\nu_k$, which is computationally expensive.

However, the incident radiation interacts with quite smooth and often completely flat dust opacities in the disk. Thus, it makes sense to ``interchange'' the order of radiative transfer and frequency integration, and to switch from a monochromatic treatment to a treatment with spectral bands.

We consider a coarse grid of frequency points $\{\nu_k~\vert~k=0,...,K\}$(e.g. $K\!=\!12$) which covers the whole SED, ranging from 100 nm to $1000~\mu$m. Instead of $B_\nu(T)$ we consider band means as

\begin{displaymath}B_k(T) = \frac{1}{\Delta\nu_k}\int_{\nu_{k-1}}^{\nu_k} B_\nu(T)~{\rm d}\nu
\end{displaymath} (29)

where $\Delta\nu_k=\nu_k-\nu_{k-1}$. In a similar way, we treat the intensities, mean intensities and opacities
$\displaystyle I_k = \frac{1}{\Delta\nu_k}\int_{\nu_{k-1}}^{\nu_k} I_\nu~{\rm d}\nu$     (30)
$\displaystyle J_k = \frac{1}{\Delta\nu_k}\int_{\nu_{k-1}}^{\nu_k} J_\nu~{\rm d}\nu$     (31)
$\displaystyle \kappa_k = \frac{1}{\Delta\nu_k}
\int_{\nu_{k-1}}^{\nu_k} \kappa_\nu~{\rm d}\nu .$     (32)

Henceforth, we exchange the index $\nu$ by the index k in all equations in this section, 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.

4.6 Dust kind, abundance, size distribution, and opacities

We assume a uniform dust abundance and size distribution throughout the disk. The dust particle density is given by

\begin{displaymath}n_{\rm d} = \int_{a_{\rm min}}^{a_{\rm max}}\!\!\! f(a)~{\rm d}a ,
\end{displaymath} (33)

where a is the particle radius and f(a) is the dust size distribution function $\rm [cm^{-4}]$, which is assumed to be given by a powerlaw as $f(a)\!=\!f_{\rm const}a^{-{a_{\rm pow}}}$. The moments of the size distribution are

\begin{displaymath}\langle a^j \rangle = \frac{1}{n_{\rm d}}
\int_{a_{\rm min}}^{a_{\rm max}}\!\!\! a^j~f(a)~{\rm d}a .
\end{displaymath} (34)

The constant in the powerlaw size distribution $f_{\rm const}$ is determined by the requirement that the dust mass density

\begin{displaymath}\rho_{\rm d} = n_{\rm d}~\rho_{\rm gr} \frac{4\pi}{3} \langle a^3 \rangle
\end{displaymath} (35)

is given by a specified fraction of the gas mas density $\rho_{\rm d}/\rho$. $\rho_{\rm gr}$ is the dust material mass density.

The dust opacities are calculated from effective medium theory (Bruggeman 1935) and Mie theory (M IEX from Wolf, according to Voshchinnikov 2002). Any uniform volume mix of solid materials with known optical constants can be used. The dust opacities are calculated as

\begin{displaymath}\kappa_\lambda^{\rm ext} = \int_{a_{\rm min}}^{a_{\rm max}}\!\!\!
\pi a^2~Q_{\rm ext}(a,\lambda)~f(a)~{\rm d}a ,
\end{displaymath} (36)

where $Q_{\rm ext}(a,\lambda)$ is the extinction efficiency. Similar formula apply for absorption and scattering opacities, $\kappa_\lambda^{\rm abs}$ and $\kappa_\lambda^{\rm sca}$, where $Q_{\rm ext}(a,\lambda)$ is replaced by $Q_{\rm abs}(a,\lambda)$ and $Q_{\rm sca}(a,\lambda)$, respectively.

Table 1:   Elements and chemical species.

5 Chemistry

The chemistry part of P ROD IM O is written in a modular form that makes it possible to consider any selection of elements and chemical species. In the models presented in this paper, we consider chemical reactions involving $N_{\rm el}\!=\!9$ elements among $N_{\rm sp}\!=\!71$ atomic, ionic, molecular and ice species as listed to Table 1.

The rate coefficients R are mostly taken from the U MIST 2006 data compilation (Woodall et al. 2007). Among the species listed in Table 1 we find 911 U MIST chemical reactions, 21 of them have multiple $T_{\hspace*{-0.2ex}\rm g}$-fits. We add 39 further reactions which are either not included in U MIST or are treated in a more sophisticated way, as explained in Sects. 5.2 to 5.5. Among the altogether 950 reactions, there are 74 photo reactions, 177 neutral-neutral and 299 ion-neutral reactions, 209 charge-exchange reactions, 46 cosmic ray and cosmic ray particle induced photo reactions, and 26 three-body reactions. The net formation rate of a chemical species i is calculated as

                           $\displaystyle \frac{{\rm d} n_i}{{\rm d}t}$ = $\displaystyle \sum\limits_{jk\ell} R_{jk\to i\ell}(T_{\hspace*{-0.2ex}\rm g})~n...
...m\limits_{j\ell} \big(R^{\rm ph}_{j\to i\ell}
+R^{\rm cr}_{j\to i\ell}\big)~n_j$  
    $\displaystyle -~n_i\!\left(\sum\limits_{jk\ell} R_{i\ell\to\!jk}~n_\ell
~+\sum\limits_{jk} \big(R^{\rm ph}_{i\to\!jk}
+R^{\rm ~cr}_{i\to\!jk}\big)\right)$ (37)

where $R_{jk\to i\ell}$ designates (two-body) gas phase reactions between two reactants j and k, forming two products i and $\ell$. $R^{\rm ph}_{i\to\!jk}$ indicates a photo-reaction which depend on the local strength of the UV radiation field, and $R^{\rm ~cr}_{j\to i\ell}$ a cosmic ray induced reaction.

5.1 Photo-reactions

Photon induced reaction rates can generally be written as

\begin{displaymath}R^{\rm ph} \;=\; 4\pi \int\!\sigma(\nu)~\frac{J_\nu}{h\nu} ~{...
...{1}{h} \int\!\sigma(\lambda)~\lambda u_\lambda ~{\rm d}\lambda
\end{displaymath} (38)

where $\sigma(\nu)$ is the photo cross-section of the reaction, $\nu$the frequency, $\lambda$ the wavelength, h the Planck constant, and $u_\lambda\!=\!\frac{4\pi}{c}J_\lambda$ the spectral photon energy density $[\rm erg~cm^{-4}]$, respectively. The proper calculation of the photo-rates $R^{\rm ph}~\rm [s^{-1}]$ according to Eq. (38) would require the calculation of a detailed (i.e. UV-line resolved) radiative transfer including molecular opacities to account for self-shielding effects (to get $J_{\nu}$) as well as detailed knowledge about the wavelength-dependent cross section  $\sigma(\lambda)$, which is not always available.

In this paper, we will apply the U MIST 2006 photo reaction rates in combination with molecular self-shielding factors from the literature instead. The application of detailed molecular UV cross sections in the calculated UV radiation field will be addressed in a future paper.

In the U MIST database, photo-rates are given as

\begin{displaymath}R^{\rm ph}_{\rm\scriptscriptstyle UMIST} \;=\; 2~
\chi_0\;\alpha~\exp(-~\gamma~A_V^{\rm\scriptscriptstyle UMIST}) ,
\end{displaymath} (39)

where $\chi_0$ is the unattenuated strength of the UV radiation field with respect to a standard interstellar radiation field, $\alpha $ the photo-rate in this standard ISM radiation field, and $A_V^{\rm\scriptscriptstyle UMIST}$ the extinction at visual wavelengths toward the UV light source. The photo-rates are derived for semi-infinite slab geometry, that is radiation is coming only from $2\pi$; this explains the factor 2 in Eq. (39). In arbitrary radiation fields, for other than 1D slab geometries, and for dust properties different from the ones used in UMIST, it is not obvious how to apply Eq. (39). In the following, we will therefore carefully explain our assumptions for the application of Eq. (39) to protoplanetary disks.

Röllig et al. (2007) relate $\chi\!=\!1$ to a ``unit Draine field'' and we will follow this idea in P ROD IM O. From the original work by Draine (1978), Draine & Bertoldi (1996) deduced

\begin{displaymath}\lambda u_\lambda^{\rm Draine} = 1.71 \cdot
...rac{31.016\lambda_3^2 - 49.913\lambda_3 + 19.897}{\lambda_3^5}
\end{displaymath} (40)

for the standard ISM UV radiation field $[\rm erg~cm^{-3}]$ where $\lambda_3\!=\!\lambda/100{\rm ~nm}$. We apply an integral definition of $\chi $ as

\begin{displaymath}\chi = \int_{91.2~{\rm nm}}^{205~{\rm nm}}
\lambda u_\lambd...
...05~{\rm nm}}
\lambda u_\lambda^{\rm Draine}~{\rm d}\lambda .
\end{displaymath} (41)

The wavelength interval boundaries have been chosen to ensure coverage of the most important photo-ionisation and photo-dissociation processes (van Dishoeck et al. 2006). Numerical integration yields $F_{\rm
Draine}\!=\!\frac{1}{h}\int_{91.2~{\rm nm}}^{205~{\rm nm}} \lambda
u_\lambda^{\rm Draine}~{\rm d}\lambda
\!=\!1.921\times10^{+8}\rm\;cm^{-2}~s^{-1}$. Adopting the wavelength boundaries 91.2 nm and 205 nm for the definition of our spectral band 1, we can directly calculate $\chi $ from our banded radiative transfer method, including scattering, by

\begin{displaymath}\chi = \frac{4\pi}{h\bar{\nu}_1} ~J_1 \Delta\nu_1 ~\Big{/}~
F_{\rm Draine} ,
\end{displaymath} (42)

where we put $\bar{\nu}_1=\sqrt{\nu_0\nu_1}$. The unshielded ISM photo-rate $\alpha $ is assumed to be given by

\begin{displaymath}\alpha \;=\; \frac{1}{2h} \int\!\sigma(\lambda)~\lambda
u_\lambda^{\rm Draine} ~{\rm d}\lambda
\end{displaymath} (43)

and the coefficient $\gamma$ in Eq. (39) can be identified as an effective, frequency-averaged opacity coefficient which contains implicit information about the frequency-range of the cross section $\sigma(\nu)$, the U MIST dust opacity and the shape of the ISM radiation field assumed. Neglecting gas extinction and assuming constant dust properties along the line of sight, the UV optical depth is

\begin{displaymath}\tau_{\rm UV} ~= \hat{\kappa}^{\rm ext}_1~N_{\rm\langle H\rangle}
\end{displaymath} (44)

where $N_{\rm\langle H\rangle}$ is the hydrogen nuclei column density toward the UV light source and $\hat{\kappa}^{\rm ext}_1$ the dust extinction coefficient per H-nucleus, averaged over spectral band 1. The ratio $A_V/\tau_{\rm
UV}$ depends on the frequency-dependence of the dust opacities as

\begin{displaymath}A_V/\tau_{\rm UV} = 2.5 \log e
\cdot\kappa^{\rm ext}_{550\rm nm}/\kappa^{\rm ext}_1 .
\end{displaymath} (45)

However, in Eq. (39) we must not use AV as calculated from our choice of dust properties! Instead, we have to use the $A_V^{\rm\scriptscriptstyle UMIST}$ depth scale as used for the compilation of the U MIST database. If we would use our AV(from dust properties in the disk), the use of $\gamma$ - containing U MIST dust properties - would internally scale it to a $\tau_{\rm UV}$ that is wrong[*]. Thus, to obtain from our UV optical depth the proper $A_V^{\rm\scriptscriptstyle UMIST}$, we need the ratio $A_V^{\rm\scriptscriptstyle UMIST}/\tau_{\rm UV}$. Since the exact U MIST-ratio is not known, we calculate it according to Eq. (45) for standard ``astronomical silicate'' grains (Draine & Lee 1984) with a size distribution $f(a)\!\propto\!a^{-3.5}$ between 0.005 $\mu$m and 0.25 $\mu$m

\begin{displaymath}A_V^{\rm\scriptscriptstyle UMIST} = 0.216\;\tau_{\rm UV} ,
\end{displaymath} (46)

whereas for larger disk dust, a value around 1 is more typical.
...ludegraphics[width=6.4cm]{1821fg06.eps} \end{tabular}\\ *[-3mm]
\end{figure} Figure 4:

Comparison of the UV radiation field strengths $\chi $ (Eq. (41)) between the simplified geometrical approach (l.h.s.), taking into account only vertical and radial dust extinction (see Eq. (47)), and the results of the full 2D radiative transfer including scattering (middle) for the $M_{\rm disk}\!=\!10^{-2}~M_\odot$ model. The two dashed contour lines show $A_V\!=\!1$ and $A_V\!=\!10$ as calculated from the minimum of the radial and vertical column densities. The r.h.s. figure shows the ratio  $\chi /\chi ^{\rm geo}$, indicating that the regions mostly affected by scattering are situated roughly between $A_V\!=\!1$ and $A_V\!=\!10$.

Open with DEXTER

Another complicated problem is how to apply Eq. (39) in disk geometry. For this purpose we introduce a geometric mean intensity as it would be present, at least approximately, if only extinction but no scattering would occur

\begin{displaymath}4\pi J_1^{\rm geo} = \Omega^\star I_1^\star {\rm e}^{-\tau_{\...
...^{\rm ISM} I_1^{\rm ISM} {\rm e}^{-\tau_{\rm UV}^{\rm
ver}} .
\end{displaymath} (47)

$I_1^\star$ and $I_1^{\rm ISM}$ are the incident band-mean stellar and interstellar intensities (see Sect. 4.3), and $\tau_{\rm UV}^{\rm rad}$ and $\tau_{\rm UV}^{\rm ver}$ are the radial (toward the star) and vertical (upwards) UV optical depth, respectively. $\Omega^\star$ is the solid angle occupied by the star and $\Omega^{\rm ISM}\!=\!4\pi-\Omega^\star$ the remainder of the full solid angle. Switching to corresponding $~\chi$ variables, we find

\begin{displaymath}\chi^{\rm geo} =
\frac{\Omega^\star}{4\pi} ~\chi_0^\star {\...
...SM}}{4\pi}~\chi_0^{\rm ISM} {\rm e}^{-\tau_{\rm UV}^{\rm ver}}
\end{displaymath} (48)

where $\chi^\star_0$ is calculated with $J_1\!=\!I_1^\star$ from Eq. (42), and $\chi_0^{\rm ISM}$ with $J_1\!=\!I_1^{\rm
ISM}$. This decomposition into two slab geometries allows us to apply Eq. (39) and calculate the photo-rates as

\begin{displaymath}R^{\rm ph} = \frac{\chi}{\chi^{\rm geo}} ~\Bigg(
...,\scriptstyle\rm ~ver}^{\rm\scriptscriptstyle UMIST}}
\end{displaymath} (49)

This approach to calculate the photo-rates according to Eq. (49) can be extended for molecular self-shielding factors (see Sect. 5.2) and states a compromise between the usual two-stream approximation and a proper treatment of UV line-resolved radiative transfer according to Eq. (38). The factor $~\chi/\chi^{\rm geo}$ corrects for the mayor shortcomings of the two-stream approximation, i.e.  the effects of scattering and the assumptions about the geometry of the radiation field made in Eq. (47). Figure 4 shows that the enhancements $~\chi/\chi^{\rm geo}$ are very close to 1 in the upper, directly irradiated layers of the disk, but may be as large as 105 at the inner rim and in the warm intermediate layer of the disk, and about 103 in the outer midplane due to scattering. $\alpha^\star$ is the unshielded photo-rate for $4\pi$-irradiation with star light, divided by $\chi^\star_0$. It can be calculated from $\sigma(\lambda)$ if known, or is assumed to be identical to $2\alpha$ otherwise.

5.2 Special UV photo reactions

For the photo-ionisation of neutral carbon, $\alpha\!=\!3\times10^{-10}~\rm s^{-1}$ is taken from the U MIST database, whereas $\alpha^\star$ is calculated according to the frequency-dependent stellar irradiation and the bound-free cross section of Osterbrock (1989). Molecular shielding by H2 and self-shielding is taken into account via the following factors from Kamp & Bertoldi (2000)
                          $\displaystyle s_{\rm C,C}$ = $\displaystyle \exp(-\sigma_{\!\rm C}^{\rm bf} N_{\rm C})$ (50)
$\displaystyle s_{\rm C,H_2}$ = $\displaystyle \exp\left(-0.9~T_{\hspace*{-0.2ex}\rm g}^{0.27}
\Big(\frac{N_{\rm H_2}}{10^{22}\rm cm^{-2}}\Big)^{0.45}\right)$ (51)

\begin{displaymath}R^{\rm ph}_{\rm C} =
s_{\rm C,C}~s_{\rm C,H_2}\;\chi_0~\alpha~\exp(-\tau_{\rm UV}) ,
\end{displaymath} (52)

i.e. we refrain from an indirect formulation with AV in cases we have the cross sections at hand. The approximation of H2 shielding for the C ionisation is strictly valid at low temperatures only ( $T\!<\!300~$K); for higher temperatures the factor 0.9 should be dropped. $N_{\rm C}$ and NH2 are the neutral carbon and molecular hydrogen column densities toward the UV light source, respectively, and $\sigma_{\!\rm C}^{\rm bf}\!=\!1.1\times10^{-17}\rm
cm^2$ the FUV-averaged cross section. The neutral carbon ionisation rates due to radial and vertical UV irradiation are calculated separately according to Eq. (52), then multiplied by the respective solid angles and added together, and then corrected for scattering as in Eq. (49).

For the photo-dissociation rate of molecular hydrogen, the same procedure applies with a H2 self-shielding factor taken from Draine & Bertoldi (1996, see their Eq. 37). We assume $\alpha~=4.2\times10^{-11}~\rm s^{-1}$ (Draine & Bertoldi 1996).

\begin{displaymath}s_{\rm H_2,H_2} =
\frac{0.965}{(1+x/b_{\rm H_2})^2}
+ \frac{0.035}{c~\exp(8.5\times 10^{-4}~c)}
\end{displaymath} (53)

with $x\!=\!N_{H_2}/5\times 10^{14}\rm ~cm^{-2}$, $b_{\rm H_2}$ the H2 UV line broadening parameter in [km/s] and $c\!=\!(1+x)^{1/2}$. The line broadening parameter b is defined as ${\it FWHM}/(4 \ln
2)^{1/2}$; it contains the sum of thermal and turbulent velocities $b
= (\frac{2 k T}{m} + \Delta v^2)^{1/2}$. Observations of line width in protoplanetary disks show that the turbulent velocities are below 0.1 km s-1 (e.g. Guilloteau & Dutrey 1998; Simon et al. 2000).

The CO photo-dissociation rate is calculated from detailed band opacities in a similar fashion, taking into account the shielding by molecular hydrogen and the self-shielding.

\begin{displaymath}s_{\rm CO,H_2} = \exp\left(-5\times 10^{\displaystyle x}\right)
\end{displaymath} (54)

with $x\!=\!9.555\times 10^{-4}\left(\log_{10}\!N_{H_2}\right)^{2.684}
- 3.976$. The CO photodissociation rate for each band is interpolated from pre-tabulated rates using the CO column density, gas temperature and line broadening parameter as input (Kamp & Bertoldi 2000; Bertoldi & Draine 1996).

5.3 H2 formation on grains

The formation of H2 on grain surfaces $\rm H + H + {\rm grain} \to
H_2 + {\rm grain}$ is taken into account according to (Cazaux & Tielens 2002)

\begin{displaymath}R_{\rm H_2} = \frac{1}{2}~v^{\rm th}_{\rm H}(T_{\hspace*{-0.2...
... a^2\rangle
~\alpha_H~\varepsilon(T_{\hspace*{-0.2ex}\rm d})
\end{displaymath} (55)

with latest updates for the temperature-dependent efficiency $\varepsilon(T_{\hspace*{-0.2ex}\rm d})$ from Cazaux (2008, priv. comm.). $v^{\rm
th}_{\rm H}\!=\!(kT_{\hspace*{-0.2ex}\rm g}/(2\pi~m_{\rm H}))^{1/2}$ is the thermal relative velocity of the hydrogen atom, $n_{\rm d}~4\pi\langle
a^2\rangle$ is the total surface of the dust component per volume, and $\alpha_H\!\approx\!0.223$ is the sticking coefficient, which results in $\frac{1}{2}~v^{\rm th}_{\rm H}(100~{\rm K})~4\pi\langle
a^2\rangle~(n_{\rm d}/n_{{\rm\langle H\rangle}})~\alpha_H\!=\!3\times 10^{-17}\rm cm^3/s$ for standard ISM grain parameter $\rho_{\rm d}/\rho\!=\!0.01$, ${a_{\rm min}}\!=\!0.005~\mu$m, ${a_{\rm max}}\!=\!0.25~\mu$m, ${a_{\rm pow}}\!=\!3.5$ and $\rho_{\rm gr}\!=\!2.5~$g/cm-3. The rate coefficient $R_{\rm
H_2}$ still needs to be multiplied by the neutral hydrogen particle density $n_{\rm H}$ to get the H2 formation rate $[\rm

5.4 Chemistry of excited H2

13 reactions for vibrationally excited molecular hydrogen H$_2^\star $are taken into account as described in Tielens & Hollenbach (1985). The FUV pumping rate ${\rm H_2} + h\nu \to {\rm H_2}^\star + h\nu'$ is assumed to be 10 times the H2 photo-dissociation rate. Two additional reactions are added for the collisional excitation by H and H2 as inverse of the de-excitation reactions

  $\textstyle {\rm H_2 + H \to H_2^\star + H} :$  $\displaystyle \quad R=C_{\rm ul}^{\rm H}(T_{\hspace*{-0.2ex}\rm g})
~\exp(-\Delta E/kT_{\hspace*{-0.2ex}\rm g})$ (56)
  $\textstyle {\rm H_2 + H_2 \to H_2^\star + H_2} :$  $\displaystyle \quad R=C_{\rm ul}^{\rm H_2}(T_{\hspace*{-0.2ex}\rm g})
~\exp(-\Delta E/kT_{\hspace*{-0.2ex}\rm g}) ,$ (57)

where the energy of the pseudo vibrational level $\Delta E=2.6~$eV as well as the collisional de-excitation rate coefficients $C_{\rm ul}^{\rm
H}$ and $C_{\rm ul}^{\rm H_2} \;\rm [cm^{-3}s^{-1}]$ are given in Tielens & Hollenbach (1985).

5.5 Ice formation and evaporation

The formation of ice mantles on dust grains plays an important role for the chemistry in the dark and shielded midplane. At the moment, five ices are considered: CO#, CO2#, H2O#, CH4# and NH3# which are treated as additional species in the chemistry (Sects. 5 and 5.6). Apart from the adsorption and desorption reactions of these species and the H2formation on grains (Sect. 5.3) no other surface reactions are currently taken into account. In particular, we do not form water on grain surfaces.

Considering collisional adsorption, and thermal, cosmic-ray and photo-desorption, the total formation rate of ice species i is

\begin{displaymath}\frac{{\rm d} n_{i\char93 }}{{\rm d}t} = n_i R^{\rm ads}_i ~-...
... des, th}_i\!
+R^{\rm des, ph}_i\!
+R^{\rm des, cr}_i\right)
\end{displaymath} (58)

where $n_{i\char93 }$ is the density of ice units i and $n_{i\char93 }^{\rm desorb}$ the fraction of $n_{i\char93 }$ located in the uppermost active surface layers of the ice mantle.

5.5.1 Adsorption

A gas species will adsorb on grain surfaces upon collision. The adsorption rate [s-1] is the product of the sticking coefficient $\alpha $, the total grain surface area per volume $4\pi\langle
a^2\rangle n_{\rm d}$ and the thermal velocity $v^{\rm
th}_i\!=\!(kT_{\hspace*{-0.2ex}\rm g}/(2\pi m_i))^{1/2}$

\begin{displaymath}R^{\rm ads}_i = 4\pi\langle a^2\rangle n_{\rm d}~\alpha~v^{\rm th}_i
\end{displaymath} (59)

where mi is the mass of gas species i. We assume unit sticking coefficient ( $\alpha\!=\!1$) for all species heavier than helium (Burke & Hollenbach 1983).

5.5.2 Desorption

A chemical species with internal energy greater than the energy that binds it to a grain surface will desorb. Desorption mechanisms depend on the source of the internal energy.

1. Thermal desorption: An ice species i at the surface of a grain at temperature $T_{\rm d}$ has probability to desorb

\begin{displaymath}R^{\rm des, th}_i = \nu^{\rm osc}_i
\exp\left(-\frac{E^{\rm ads}_i}{kT_{\rm d}}\right) ,
\end{displaymath} (60)

where $\nu^{\rm osc}_i\!=\!(2~n_{\rm surf}~k E^{\rm ads}_i/(\pi^2
m_i))^{1/2}$ is the vibrational frequency of the species in the surface potential well of ice species i, $n_{\rm surf}\!=\!1.5\times
10^{15}$ cm-2 is the surface density of adsorption sites and $E^{\rm ads}_i$ is the adsorption binding energy. The adopted values are provided in Table 2. Following Aikawa et al. (1996), the number density of ice units at the active surface $n_{i\char93 }^{\rm desorb}$ is given by

\begin{displaymath}n_{i\char93 }^{\rm desorb} = \left(\begin{array}{ll}
...d\ \ n^{\rm ice}_{\rm tot}\ge n_{\rm act}
\end{displaymath} (61)

where $n_{\rm act}\!=\!4\pi\langle a^2\rangle n_{\rm d}~n_{\rm surf}N_{\rm Lay}$is the number of active surface places in the ice mantle per volume and $N_{\rm Lay}$ is the number of surface layers to be considered as ``active''. We assume $N_{\rm Lay}\!=\!2$ in accordance with Aikawa et al. (1996). $n^{\rm ice}_{\rm tot}\!=\!\sum_j n_{j\char93 }$ is the number density of ices.

2. Photo-desorption: Absorption of a UV photon by a surface species can increase the species internal energy enough to induce desorption. The photo-desorption rate of species i is given by

\begin{displaymath}R^{\rm des, ph}_i = \pi\langle a^2\rangle \frac{n_{\rm d}}{n_{\rm act}}
~Y_i~\chi F_{\rm Draine}
\end{displaymath} (62)

where Yi is the photo-desorption yield (see Table 2), $\chi F_{\rm Draine}$ is a flux-like measure of the local UV energy density [photons/cm2/s] computed from continuum radiative transfer (Eqs. (41), (42)). Photo-desorption can enhance gas-phase water abundances by orders of magnitude in outer region of disks (Willacy & Langer 2000; Öberg et al. 2008; Dominik et al. 2005).

Table 2:   Adsorption energies and photo-desorption yields.

3. Cosmic-ray induced desorption: Cosmic-rays hitting a grain can locally heat the surface and trigger desorption. Cosmic-rays can penetrate deep into obscured regions, maintaining a minimum amount of species in the gas-phase. Cosmic-ray fluxes in disks may be higher than in molecular clouds because of the stellar energetic particles in addition to the galactic component. X-ray photons can also penetrate deep inside the disk and locally heat a dust grain but X-ray induced desorption is not included in the code yet. We adopt for the cosmic-ray desorption the formalism of Hasegawa & Herbst (1993).

\begin{displaymath}R^{\rm des, cr}_i = f(70~{\rm K})~R^{\rm des, th}_i(70~{\rm K})~
\frac{\zeta_{\rm CR}}{5\times 10^{-17}~{\rm s}^{-1}}
\end{displaymath} (63)

where $\zeta_{\rm CR}$ is the cosmic ray ionisation rate of H2, $f(70~{\rm K})\!=\!3.16 \times 10^{-19}$ the ``duty-cycle'' of the grain at 70 K and $R^{\rm des, th}_i(70~\rm K)$ the thermal desorption rate for species i at temperature $T_{\hspace*{-0.2ex}\rm d}\!=\!70$ K. The adopted value for $f(70~\rm K)$ is strictly valid only for 0.1 $\mu$m grains in dense molecular clouds.

5.6 Kinetic chemical equilibrium

Assuming kinetic chemical equilibrium in the gas phase, and between gas and ice species, we have $\frac{{\rm d} n_i}{{\rm d}t}\!=\!0$ in Eq. (37) and obtain $i\!=\!1~...~N_{\rm sp}$ non-linear equations for the unknown particle densities nj $(j\!=\!1~...~N_{\rm sp})$

\begin{displaymath}F_{\!i}(n_j) ~=~0 .
\end{displaymath} (64)

It is noteworthy that the electron density $n_{\rm e}$ is not among the unknowns, but is replaced by the constraint of charge conservation

\begin{displaymath}n_{\rm e} = \sum_j n_j ~z_j
\end{displaymath} (65)

where zj is the charge of species j in units of the elementary charge. The explicit dependency of $n_{\rm e}$ on the particle densities nj causes additional entries in the chemical Jacobian d $ F_i/{\rm d} n_j = \partial
F_j/\partial n_j + \partial F_i/\partial n_{\rm e} \cdot \partial n_{\rm
e}/\partial n_j$.

5.7 Element conservation

The system of Eqs. (37) is degenerate because every individual chemical reaction obeys several element conservation constraints, and therefore, certain linear combinations of $F_{\!j}$can be found which cancel out, making the equation system under-determined. Only if the element conservation is implemented in addition, the system (Eqs. (37)) becomes well-defined.

Considering the total hydrogen nuclei density $n_{{\rm\langle H\rangle}}$ as given, the conservation of element k can be written as

\begin{displaymath}n_{{\rm\langle H\rangle}}~\epsilon_k - \sum_i n_i~\nu_{i,k} ~=~ 0 ,
\end{displaymath} (66)

resulting in $N_{\rm el}$ auxiliary conditions. $\epsilon_k$ is the elemental abundance of element k normalised to hydrogen and $\nu_{i,k}$ are the stoichiometric coefficient of species i with respect to element k.

Alternatively, the gas pressure p may be considered as the given quantity and the relative element conservation can be expressed by

\begin{displaymath}\hat\epsilon_k \Big(n_{\rm e}~m_{\rm e}+\sum_i n_i m_i\Big)
- \sum_i n_i~\nu_{i,k}~m_k ~=~0 ,
\end{displaymath} (67)

where $\hat\epsilon_k=\epsilon_k m_k/(\sum_{k'} \epsilon_{k'} m_{k'})$is the relative mass fraction of element k, mi the mass of a gas particle of kind i and mk the mass of element k. Since summing up all Eqs. (67) for $k\!=\!1~...~N_{\rm el}$ results in $\rho\!=\!\rho$, one of these equations is redundant and can be replaced by the constraint of given pressure

\begin{displaymath}p - \sum_i n_i \left( 1 + z_i\right) kT_{\hspace*{-0.2ex}\rm g}~=~0 ,
\end{displaymath} (68)

where k the Boltzmann constant. The element conservation is implemented by replacing  $N_{\rm el}$ selected components of $F_{\!j}$in Eq. (64) by these auxiliary conditions, either according to Eq. (66) or according to Eqs. (67) and (68), after suitable normalisation. For this purpose, we choose for every element k the index j that belongs to the most abundant species containing this element. Equation (68) overwrites the entry for the most abundant H-containing species.

The global iteration, which solves the hydrostatic disk structure consistently with the chemistry and heating & cooling balance (see Fig. 1), is found to converge only if the chemistry is solved for constant pressure p. Since the vertical hydrostatic condition (Eq. (2)) is a pressure constraint, it is essential to ensure that the chemistry solver, coupled to the $T_{\hspace*{-0.2ex}\rm g}$-determination via heating ${\rm\hspace*{0.7ex}\&\hspace*{0.7ex}}$cooling balance, is not allowed to change p as it would be the case if  $n_{{\rm\langle H\rangle}}$ was fixed. At given pressure p, $T_{\hspace*{-0.2ex}\rm g}$ may be found to increase during the course of the iteration, but only if simultaneously  $n_{{\rm\langle H\rangle}}$ drops, thereby conserving the p-structure within one global iteration step.

5.8 Numerical solution of chemistry

The non-linear equation system (64), expressing the kinetic chemical equilibrium including element conservation, is usually solved by means of a self-developed, globally convergent Newton-Raphson method. A quick and reliable numerical solution of Eqs. (64) is crucial for the computational time consumption, stability, and global convergence of our model. Our numerical experience shows that a careful storage of converged results (particle densities) is the key to increase stability and performance. These particle densities are used as initial guesses for the next time the Newton-Raphson method is invoked, either in form of a downward-outward sweep through the grid (first iteration), or from the last results of the same point (following iterations).

In cases, where the solution by the Newton-Raphson method fails, we fall back to the time-dependent case and solve Eqs. (37) by means of the ODE solver L IMEX (Deuflhard & Nowak 1987) for 107 yr, which is much slower but in practice gives the same results as the Newton-Raphson method.

6 Gas thermal balance

The net gain of thermal kinetic energy is written as

\begin{displaymath}\frac{de}{dt} =
\sum\limits_k \Gamma_k(T_{\hspace*{-0.2ex}\...
... \sum\limits_k \Lambda_k(T_{\hspace*{-0.2ex}\rm g},n_{\rm sp})
\end{displaymath} (69)

where $\Gamma_k$ and $\Lambda_k$ are the various heating and cooling rates $\rm [erg~cm^{-3}~s^{-1}]$ which are detailed in the forthcoming sub-sections. Restricting ourselves to the case of thermal energy balance, we assume d $e/{\rm d}t\!=\!0$ in the following and Eq. (69) states an implicit equation for the unknown kinetic gas temperature $T_{\hspace*{-0.2ex}\rm g}$. Since the heating and cooling rates depend not only on $T_{\hspace*{-0.2ex}\rm g}$, but also on the particle densities $n_{\rm
sp}$, which themselves depend on $T_{\hspace*{-0.2ex}\rm g}$, an iterative process is required during which $T_{\hspace*{-0.2ex}\rm g}$ is varied and the the chemistry is re-solved until $T_{\hspace*{-0.2ex}\rm g}$ satisfies Eq. (69).

6.1 Non-LTE treatment of atoms, ions and molecules

The most basic interaction between matter and radiation is the absorption and emission of line photons by a gas particle, which can be an atom, ion or molecule. We consider a N-level system with bound-bound transitions only and calculate the level populations $n_j~\rm [cm^{-3}]$ by means of the statistical equations

\begin{displaymath}n_i \sum\limits_{j\ne i} R_{ij} = \sum\limits_{j\ne i} n_j R_{ji}
\end{displaymath} (70)

which are solved together with the equation for the conservation of the total particle density of the considered species $\sum_i n_i\!=\!n_{\rm sp}$. The rate coefficients are given by (Mihalas 1978):
$\displaystyle R_{\rm ul} = A_{\rm ul}
+ B_{\rm ul}\overline{J}_{\rm ul} + C_{\rm ul}$      
$\displaystyle R_{lu} = B_{lu}\overline{J}_{\rm ul} + C_{lu},$     (71)

where u and l label an upper and lower level, respectively. $A_{\rm ul}$, $B_{\rm ul}$, Blu, $C_{\rm ul}$ and Clu are the Einstein coefficients for spontaneous emission, absorption, stimulated emission and the rate coefficients for collisional (de-)excitation, respectively. Additionally we have the Einstein relations $B_{\rm ul}/A_{\rm ul}\!=\!c^2/(2h\nu_{\rm ul}^3)$, $B_{lu}/B_{\rm ul}\!=\!g_u/g_l$and the detailed balance relation $C_{lu}/C_{\rm ul}\!=\!g_u/g_l \cdot
\exp(-\Delta E_{\rm ul}/kT_{\hspace*{-0.2ex}\rm g})$, where $\nu_{\rm ul}$, gu, gl and $\Delta
E_{\rm ul}$ are the line centre frequency, the statistical weights of the upper and lower level and the energy difference, respectively. The line integrated mean intensity is given by

\begin{displaymath}\overline{J}_{\rm ul} = \frac{1}{4\pi}
\iint\!\phi_{\rm ul}(\nu,\vec{n})~I_{\nu}(\vec{n})~{\rm d}\nu~\rm d\Omega
\end{displaymath} (72)

where $\phi_{\rm ul}(\nu,\vec{n})$ is the line profile function in direction $\vec{n}$.

6.1.1 Escape-probability treatment

\end{figure} Figure 5:

Different pumping and escape probabilities according to the predominantly radial irradiation and the predominantly vertical escape.

Open with DEXTER

The spectral intensity $I_\nu$ in Eq. (72) is affected by line absorption and emission. Assuming that the line source function (Eq. (74)) varies slowly in a local environment where the line optical depths (Eq. (75)) build up rapidly, we can approximate for a static, plane-parallel medium

\begin{displaymath}I_{\nu}(\mu) ~\approx~
I^{\rm cont}_{\nu_{\rm ul}}(\mu)\exp...
...c{\phi_{\rm ul}(\nu)~\tau^{\rm ver}_{\rm ul}}{\mu}\Big)\right)
\end{displaymath} (73)

where $I^{\rm cont}_{\nu_{\rm ul}}(\mu)$ is the continuous background intensity which propagates backward along the ray in direction $\mu$. The direction $\mu=1$ points ``outward'' ( $\mu\!=\!\cos\theta$). The line source function and the perpendicular line optical depth are given by (Mihalas 1978)
                          $\displaystyle S^{\rm L}_{\rm ul}$ = $\displaystyle \frac{2h\nu_{\rm ul}^3}{c^2}
\left( \frac{g_u n_l}{g_l n_u} - 1 \right)^{-1}$ (74)
$\displaystyle \tau^{\rm ver}_{\rm ul}$ = $\displaystyle \frac{A_{\rm ul}~c^3}{8\pi\nu_{\rm ul}^3\Delta {\rm v}_D}
\int_z^{z_{\rm max}}\!\!\Big(n_l(z')\frac{g_u}{g_l}-n_u(z')\Big)~{\rm d}z'$ (75)

where $\Delta {\rm v}_D$ is the (turbulent + thermal) velocity Doppler width of the line, assumed to be constant along the line of sight in Eq. (75). Equations (72) and (73) can be combined to find
                                   $\displaystyle \overline{J}_{\rm ul}$ = $\displaystyle \frac{1}{2} \int_{-1}^1
p^{\rm ~esc}_{\rm ul}(\mu)~I^{\rm cont}_{...
+ \left(1-p^{\rm ~esc}_{\rm ul}(\mu)\right) S^{\rm L}_{\rm ul} \;{\rm d}\mu$ (76)
  $\textstyle \approx$ $\displaystyle P^{\rm ~pump}_{\rm ul}~J^{\rm cont}_{\nu_{\rm ul}}
\;+\; \left(1-P^{\rm ~esc}_{\rm ul}\right)~S^{\rm L}_{\rm ul}$ (77)

where the direction-dependent and the mean escape probabilities are found to be
                          v$\displaystyle p^{\rm ~esc}_{\rm ul}(\mu)$ = $\displaystyle \int_{-\infty}^{+\infty} \phi(x)~
\exp\big(-\frac{\tau^{\rm ver}_{\rm ul}~\phi(x)}{\mu}\big)~{\rm d}x$ (78)
$\displaystyle P^{\rm ~esc}_{\rm ul}$ = $\displaystyle \frac{1}{2}\int_{-1}^{~1} p^{\rm ~esc}_{\rm ul}(\mu)\;{\rm d}\mu$ (79)

with dimensionless line profile function $\phi(x)\!=\!\exp(-x^2)/\!\sqrt{\pi}$, $x\!=\!(\nu-\nu_{\rm ul})/\Delta\nu_D$ and frequency width $\Delta\nu_D\!=\!\nu_{\rm ul}\Delta {\rm v}_D/c$. Using Eq. (77), it is straightforward to show that the unknown line source function $S^{\rm L}_{\rm ul}$ can be eliminated, and the leading term $n_u A_{\rm ul}$ cancels out, when considering the net rate $n_u A_{\rm ul} + (n_u B_{\rm ul} - n_l B_{lu})~\overline{J}_{\rm ul} = n_u
... (n_u B_{\rm ul} - n_l
B_{lu})P^{\rm ~pump}_{\rm ul}J^{\rm cont}_{\nu_{\rm ul}}$. Thus, we can solve the statistical rate Eqs. (70) with modified rate coefficients
                          $\displaystyle \tilde{R}_{\rm ul}$ = $\displaystyle A_{\rm ul}P^{\rm ~esc}_{\rm ul}
+ B_{\rm ul}P^{\rm ~pump}_{\!ul}J^{\rm cont}_{\nu_{\rm ul}} + C_{\rm ul}$  
$\displaystyle \tilde{R}_{lu}$ = $\displaystyle B_{lu}P^{\rm ~pump}_{\!ul}J^{\rm cont}_{\nu_{\rm ul}} + C_{lu}
,$ (80)

which is known as escape-probability formalism (Avrett & Hummer 1965; Mihalas 1978). $P^{\rm ~esc}_{\!ul}$ is the mean probability for line photons emitted from the current position to escape the local environment and $P^{\rm ~pump}_{\!ul}$ the mean probability for continuum photons to arrive at the current position. $J^{\rm cont}_{\nu_{\rm ul}}$ is the continuous mean intensity at line centre frequency $\nu_{\rm ul}$, as would be present if no line transfer effects took place. In semi-infinite slab symmetry, all directions $\mu\!<\!0$ have infinite line optical depth and can be discarded from the calculation of the escape probabilities
                            $\displaystyle P^{\rm ~esc}_{\rm ul}(\tau^{\rm ver}_{\rm ul})$ = $\displaystyle \frac{1}{2}\int_{-\infty}^{+\infty}\!\!\!\!
\exp\big(-\frac{\tau^{\rm ver}_{\rm ul}\phi(x)}{\mu}\big)\;{\rm d}\mu~{\rm d}x$ (81)
  = $\displaystyle \frac{1}{2}\int_{-\infty}^{+\infty}\!\!\!\!
\phi(x)~E_2\big(\tau^{\rm ver}_{\rm ul}\phi(x)\big)~{\rm d}x.$ (82)

This function is numerically fitted as

\begin{displaymath}P^{\rm ~esc}_{\rm ul}(\tau) = \left\{ \begin{array}{ll}
0.5 ...
...(0.4799\tau)\big)^{-0.4195} , & 9\!<\!\tau.

Considering the pumping probability as defined by Eq. (77), it is noteworthy that $P^{\rm ~pump}_{\rm ul}\!\approx\!P^{\rm ~esc}_{\rm ul}$ is only valid in an almost isotropic background radiation field. In disk symmetry, much of the pumping is due to direct star light (see Fig. 4) which has a very pointed character. In the optically thick midplane, the continuous radiation field is almost isotropic, but here the pumping is pointless, because the radiation is thermalised and the collisional processes dominate. Considering near to far IR wavelengths at a certain height above the midplane, the irradiation from underneath plays a role, but these directions are just the opposite of what is considered in Eq. (82), and so using $P^{\rm ~pump}_{\rm ul}\!\approx\!P^{\rm ~esc}_{\rm ul}$ would be strongly misleading. Thus, we approximate

$\displaystyle P^{\rm ~pump}_{\rm ul}(\tau^{\rm rad}_{\rm ul}) = \int_{-\infty}^{+\infty}\!\!\!\!
\phi(x)~\exp\big(-\tau^{\rm rad}_{\rm ul}\phi(x)\big)\;{\rm d}x$     (83)

with $\tau^{\rm rad}_{\rm ul}$ now being the radially inward line optical depth. This function is numerically fitted as

\begin{displaymath}P^{\rm ~pump}_{\rm ul}(\tau) = \left\{ \begin{array}{ll}
1 ,...
...n(0.3367\tau)\big)^{-0.4306} ,& 9\!<\!\tau.

\end{figure} Figure 6:

Continuum mean intensities as input for non-LTE modelling. The calculated band-mean mean intensities are shown for one particular point (r,z) 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. The red line shows the band-mean incident stellar intensities $\nu I_\nu ^\star $ and the blue line shows the incident interstellar intensities $\nu I_\nu ^{\rm ISM}$. The radiation field has two major components, the dust attenuated UV - near IR part, originating mainly from the star, and the self-generated mid - far IR part, originating from thermal dust emission in the disk.

Open with DEXTER

6.1.2 Background radiation field

The continuum background mean intensities $J^{\rm cont}_{\nu_{\rm ul}}$ have an important impact on the gas energy balance. For example, in strong continuum radiation fields, the reverse process to line cooling, namely line absorption followed by collisional de-excitation, dominates. $J^{\rm cont}_{\nu_{\rm ul}}$ is identified to be just given by the mean intensities calculated from the dust continuum radiative transfer (see Sect. 4). In order to obtain the required monochromatic mean continuum intensities at the line centre positions, we apply a cubic spline interpolation to the calculated local continuum $J_\nu^{\rm cont}(r,z)$ in frequency space as depicted in Fig. 6.

6.1.3 Solving the statistical equations

Equations ((70), (75), (80)) form a system of coupled equations for the unknown population numbers ni. Since the line optical depths (Eq. (75)) depend partly on the local ni, these equations must be solved iteratively. We apply a fully implicit integration scheme for the numerical solution of (Eq. (75)) where the line optical depth increment along the last downward integration step, between the previous and the current grid point, is given by the local populations, i.e. 

\begin{displaymath}\Delta\tau^{\rm ver}_{i} = \frac{A_{\rm ul}~c^3}{8\pi\nu_{\rm...
...\langle H\rangle}(z_i)-N_{\rm\langle H\rangle}(z_{i+1})\big) .
\end{displaymath} (84)

where $N_{\rm\langle H\rangle}(z_i)$ is the vertical hydrogen column density at grid point zi, and nu, nl and $n_{{\rm\langle H\rangle}}$ refer to the current grid point i. A simple $\Lambda$-type iteration scheme is found to converge within typically 0 to 20 iterations. The outward radial line optical depth increments $\Delta\tau_j^{\rm rad}$ between rj-1 and rj are calculated in a similar fully implicit fashion.

6.1.4 Calculation of the heating/cooling rate

Once the statistical equations (Eqs. (70)) have been solved, the radiative heating and cooling rates can be determined. There are two valid approaches. For the net cooling rate, one can either calculate the net creation rate of photon energy (radiative approach), or one can calculate the total destruction rate of thermal energy (collisional approach).

                           $\displaystyle \Gamma_{\rm rad}$ = $\displaystyle \sum\limits_{u>l} n_l \Delta E_{\rm ul}
P^{\rm ~pump}_{\rm ul}B_{lu}J^{\rm cont}_{\nu_{\rm ul}}$ (85)
$\displaystyle \Lambda_{\rm rad}$ = $\displaystyle \sum\limits_{u>l} n_u \Delta E_{\rm ul}
~\big(P^{\rm ~esc}_{\rm ul}A_{\rm ul} + P^{\rm ~pump}_{\rm ul}B_{\rm ul}J^{\rm cont}_{\nu_{\rm ul}}\big)$ (86)
$\displaystyle \Gamma_{\rm col}$ = $\displaystyle \sum\limits_{u>l} n_u C_{\rm ul} \Delta E_{\rm ul}$ (87)
$\displaystyle \Lambda_{\rm col}$ = $\displaystyle \sum\limits_{u>l} n_l C_{lu} \Delta E_{\rm ul}.$ (88)

Both approaches must yield the same net result $\Gamma_{\rm
rad}\!-\!\Lambda_{\rm rad}\!=\!\Gamma_{\rm coll}\!-\!\Lambda_{\rm
coll}$, which can be used to check the quality of the numerical solution. In practise, one pair of these heating/cooling rates is often huge in comparison to the other pair, e.g.  $\{\Gamma_{\rm
rad},\Lambda_{\rm rad}\} \gg \{\Gamma_{\rm col},\Lambda_{\rm col}\}$in a thin gas with radiatively controlled populations, and $\{\Gamma_{\rm rad},\Lambda_{\rm rad}\} \ll \{\Gamma_{\rm
col},\Lambda_{\rm col}\}$ in a dense gas with close to LTE populations. Thus, it is numerically favourable to choose
$\displaystyle \Gamma = \left\{\begin{array}{lccl}
\Gamma_{\rm col}, &&& \Gamma_...
\Gamma_{\rm rad}, &&& \Gamma_{\rm rad}\le \Gamma_{\rm col}
\end{array}\right.$     (89)
$\displaystyle \Lambda = \left\{\begin{array}{lccl}
\Lambda_{\rm col},&&& \Gamma...
\Lambda_{\rm rad},&&& \Gamma_{\rm rad}\le \Gamma_{\rm col}.
\end{array}\right.$     (90)

Table 3:   Assumed element abundances in (gas + ice)

Table 4:   Non-LTE model atoms, ions and molecules.

6.1.5 Atomic and molecular data

The atomic and molecular data for O I, C I, C II and H2O (energy levels Ei, statistical weights gi, Einstein coefficients $A_{\rm ul}$, and collision rates $C_{\rm ul}$ are taken from Leiden's LAMBDA-database (Schöier et al. 2005), see Table 4. In addition to these low-temperature coolants, we have included several ions as high-temperature coolants from the CHIANTI-database (Dere et al. 1997): Mg  II, Fe  II, Si  II and S  II, taking into account all energy levels up to about 60 000 cm-1. This database has collisional data for free electrons only, but since we consider only ions of abundant elements here, the electron concentration is always rather high wherever these ions are abundant. Since the electron collisional rates are typically 104 times larger than those of heavy particles, the thereby introduced error seems acceptable.

For CO, we have merged level and radiative data (Ei, gi and $A_{\rm ul}$) of the rotational and ro-vibrational states ( $v\!=\!0,1,2,3,4$) from the H ITRAN database (Rothman et al. 2005) with collisional data among the rotational levels from the LAMBDA database. For the vibrational collisions we use the $C_{1\to0}$ data for H and H2 de-exciting collisions from Neufeld & Hollenbach (1994) and for He collisions from Millikan & White (1964). The de-exciting rate coefficients for other than $1\!\to\!0$vibrational transitions are estimated according to the formula provided by Elitzur (1983)

\begin{displaymath}C_{v'\to v} = (v-v')~C_{1\to0} \exp\left(-\frac{(v-v'-1)~1.5~...
....2ex}\rm g}}
{1+1.5~\theta/T_{\hspace*{-0.2ex}\rm g}}\right)
\end{displaymath} (91)

where $\theta\!=\!\hbar\omega/k$ is the difference between the first vibrationally excited and the ground state energy. For detailed ro-vibrational modelling, these total vibrational collisional rates still need to be spread over the rotational sub-states. For simplicity, however, we assume $C_{v'J'\to vJ}\approx C_{v'\to
v}$ for every rotational state J.

For H2, the level and radiative data (quadrupole transitions) are taken from Wolniewicz et al. (1998). We include calculated collisional excitations by H (Wrathmall et al. 2007), ortho- and para-H2, and Helium (Le Bourlot et al. 1999). The H2 and H2O ortho to para abundance ratios are assumed to be at thermal equilibrium according to the gas temperature.

6.2 Specific heating processes

Below, we list further heating processes that are not covered by Sect. 6.1. Photoelectric heating, cosmic ray ionisation, carbon photo-ionisation and H2 photo-dissociation are still radiative processes, while other heating mechanisms are of chemical nature, such as H2formation heating, or of dynamical nature, such as viscous heating.

6.2.1 Photo-electric heating

UV photons impinging on dust grains can eject electrons with super-thermal velocities which then thermalise through collisions with the gas. The efficiency of this process decreases strongly with grain charge (positively charged grains are less efficient heaters). The grain charge is set by the balance of incoming UV flux that ejects electrons and collisional recombination. The collision rate for recombination scales with electron density, thermal velocity and the ratio between potential and thermal energy ( $\Phi\!=\!eU/kT_{\hspace*{-0.2ex}\rm g}$, with U being the grain potential). Thus the grain charge can be parameterised by a 'so-called' grain charge parameter (Bakes & Tielens 1994)

\begin{displaymath}x = \sqrt{T_{\hspace*{-0.2ex}\rm g}}~\frac{\chi}{n_{\rm e}} \cdot
\end{displaymath} (92)

The probability of electron ejection after photon absorption (yield), is generally taken from experimental data on bulk material with large flat surfaces, and then applied to (smaller) astrophysical dust grains to compute the photoelectric heating rates. The heating process is thought to be less effective for micron-sized grains compared to small ISM dust grains. The reason is that a photo-electron can more easily be trapped within the matrix of a large grain, thus lowering the photoelectric yield. Experimental data on realistic astrophysical dust grains is sparse and only recently (Abbas et al. 2006) carried out experiments with sub-micron to micron sized individual dust grains. They measure yields that are larger than those of bulk flat surfaces and they find an increasing yield with increasing grain size. However, the underlying physics of these experiments are not yet well understood.

Kamp & Bertoldi (2000) provide a formula to approximate the photoelectric heating rate for large graphite and silicate grains using the photoelectric yields of bulk material from Feuerbacher & Fitton (1972). For silicate grains, the photoelectric heating rate $\Gamma_{\rm PE}$and the efficiency $\epsilon$ are

                         $\displaystyle \Gamma_{\rm PE}$ = $\displaystyle 2.5\times 10^{-4} ~\sigma^{\rm abs}n_{{\rm\langle H\rangle}}~\epsilon~\chi$ (93)
$\displaystyle \epsilon$ = $\displaystyle \frac{0.06}{1+1.8\times10^{-3}~x^{~0.91}}
+\frac{y~(10^{-4} T_{\hspace*{-0.2ex}\rm g})^{1.2}}{1+0.01 x}$ (94)
y = $\displaystyle \left\{ \begin{array}{ll}
0.7 \hspace*{3mm} & , x \leq 10^{-4}\\  [1mm]
0.36 & , 10^{-4} < x \leq 1\\  [1mm]
0.15 & , x > 1
\end{array}\right.$ (95)

valid for electron particle densities $10^{-5}{\rm cm^{-3}}\!<\!n_{\rm
e}\!<\!10^5{\rm cm^{-3}}$, gas temperatures $10~{\rm
K}\!<\!T_{\hspace*{-0.2ex}\rm g}\!<10000~$K, and strength of FUV radiation field $10^{-5}\!<\!\chi\!<\!10^5$. Here, $\sigma^{\rm abs}$ is the grain absorption cross section per H-nucleus ( $\sigma^{\rm abs}n_{{\rm\langle H\rangle}}\!=\!\kappa_1^{\rm abs}$, see Eqs. (32) and (36)).

6.2.2 PAH heating

Very small dust grains such as polycyclic aromatic hydrocarbons (PAHs) are an extremely efficient heating source for the gas. The photoelectric heating rate can be written separately from the rest of the grain size distribution using only the first term of the (Bakes & Tielens 1994) efficiency formulation

\begin{displaymath}\Gamma_{\rm PE} = f_{\rm PAH}~10^{-24}n_{{\rm\langle H\rangle}}~\epsilon~\chi
\end{displaymath} (96)

where the efficiency $\epsilon$ is given by

\begin{displaymath}\epsilon = \frac{0.0487}{1 + 4\times 10^{-3}x^{~0.73}}\cdot
\end{displaymath} (97)

In the ISM, the abundance of PAHs is $f_{\rm PAH}\!=\!1$. For disks, this value can be scaled according to the observed strength of the PAH bands.

6.2.3 Carbon photo-ionisation

Ionisation of carbon releases electrons with energies around 1 eV (Black 1987). Subsequent collisions heat the gas as

\begin{displaymath}\Gamma_{\rm C} = 1.602\times 10^{-12}~R^{\rm ph}_{\rm C}~n_{\rm C}
\end{displaymath} (98)

where the photo-ionisation rate $R^{\rm ph}_{\rm C}$ is given by Eq. (52).

6.2.4 H2 photo-dissociation heating

Photo-dissociation of molecular hydrogen occurs via UV line absorption into an electronically excited state followed by spontaneous decay into an unbound state of the two hydrogen atoms. The kinetic energy of such H-atoms is typically 0.4 eV (Stephens & Dalgarno 1973), leading to an approximate heating rate of

\begin{displaymath}\Gamma_{\rm ph}^{\rm H_2} = 6.4 \times 10^{-13}~R_{\rm ph}^{\rm H_2}
~n_{\rm H_2}.
\end{displaymath} (99)

Here, $R_{\rm ph}^{\rm H_2}$ is the H2 photo-dissociation rate given in Sect. 5.2 including dust and H2self-shielding.

6.2.5 cosmic ray heating

Cosmic rays have a typical attenuation depth of 96 g cm-2 and thus reach much deeper than stellar FUV photons ($\sim$10-3 g cm-2, see Bergin et al. 2007, for an overview). They ionise atomic and molecular hydrogen and this inputs approximately 3.5 and 8 eV into the gas for H and H2, respectively (Jonkheid et al. 2004). The heating rate can then be written as

\begin{displaymath}\Gamma_{\rm CR} = \zeta_{\rm CR} \left( 5.5 \times 10^{-12} n_{\rm H}
+ 2.5 \times 10^{-11} n_{\rm H_2} \right)
\end{displaymath} (100)

where $\zeta_{\rm CR}$ is the primary cosmic ray ionisation rate.

6.2.6 H2 formation heating

The formation of H2 on dust surfaces releases the binding energy of 4.48 eV. Due to the lack of laboratory data, we follow the approach by Black & Dalgarno (1976) and assume that this energy is equally distributed over translation, vibration and rotation. Hence, about 1.5 eV per reaction is liberated as heat

\begin{displaymath}\Gamma_{\rm form}^{\rm H_2} = 2.39 \times 10^{-12}~R_{\rm H_2}~n_{\rm H}
\end{displaymath} (101)

where the H2 formation rate $R_{\rm
H_2}$ is given in Sect. 5.3.

6.2.7 Heating by collisional de-excitation of H$_2^\star $

The fluorescent excitation of H2 by UV photons $\rm H_2 + h\nu \to
H_2^{\star\star} \to \rm H_2^\star + h\nu'$ produces vibrationally excited molecular hydrogen H$_2^\star $ (Tielens & Hollenbach 1985), and the vibrational excitation energy can be converted into thermal energy by collisions. The heating rate is

                             $\displaystyle \Gamma_{\rm coll}^{\rm H_2^\star}$ = $\displaystyle \Delta E~
R^{\rm coll}_{\rm H_2\to H_2^\star}
\left( n_{\rm H_2^\...
...- n_{\rm H_2}
\exp\Big(-\frac{\Delta E}{kT_{\hspace*{-0.2ex}\rm g}}\Big)\right)$ (102)
$\displaystyle R^{\rm coll}_{\rm H_2\to H_2^\star}$ = $\displaystyle n_{\rm H} ~C_{\rm ul}^{\rm H}(T_{\hspace*{-0.2ex}\rm g})
+ n_{\rm H_2}~C_{\rm ul}^{\rm H_2}(T_{\hspace*{-0.2ex}\rm g})$ (103)

where the excitation energy of the pseudo vibration level $\Delta E$ and the collisional de-excitation rates $C_{\rm ul}$ are given in Tielens & Hollenbach (1985), see also Sect. 5.4. The second term in Eq. (102) corrects for collisional excitation.

6.2.8 Viscous heating

Due to high optical thickness, radiative heating cannot penetrate efficiently to the midplane. These dense layers can instead also be heated by local viscous dissipation (Frank et al. 1992)

\begin{displaymath}\Gamma_{\rm vis} = \frac{9}{4}~\rho~\nu_{\rm kin}~\Omega_{\rm kep}^2
\end{displaymath} (104)

In the absence of a well-understood mechanism, angular momentum transport is conceptualised using the kinematic - or turbulent - viscosity $\nu_{\rm kin}$ often parameterised as an $\alpha $-viscosity (Shakura & Syunyaev 1973)

\begin{displaymath}\nu_{\rm kin} = \alpha~c_T~H_{\rm g} ,
\end{displaymath} (105)

where $\alpha $ is a dimensionless scaling factor, $c_T^2\!=\!p/\rho$is the isothermal sound speed, $H_{\rm g}\!=\!c_T/\Omega_{\rm kep}$ is the gas scale height, and $\Omega_{\rm kep}\!=\!v_\phi/r$ is the Keplerian angular velocity (see Eq. (4)). For $r\!\la\!2$ AU, the viscous heating is known to be capable of dominating the energy balance in the midplane (D'Alessio et al. 1998)[*].

6.3 Specific cooling processes

Most cooling processes are radiative in nature and covered in Sect. 6.1. However, two prominent high temperature cooling processes are treated in a simpler approximative fashion: Lyman-$\alpha $ and O I-630 nm cooling.

6.3.1 Ly-$\alpha $ cooling

Cooling through the Lyman-$\alpha $ line becomes efficient at temperatures of a few 1000 K (Sternberg & Dalgarno 1989). Given the densities of atomic hydrogen $n_{\rm H}$ and electrons $n_{\rm e}$, the cooling rate can be written as

\begin{displaymath}\Lambda_{\rm Ly-\alpha} = 7.3 \times 10^{-19} n_{\rm H}~n_{\rm e}
\exp\left( -118~400/T_{\hspace*{-0.2ex}\rm g}\right) .
\end{displaymath} (106)

6.3.2 O I-630nm cooling

Line emission from the meta-stable 1D level of neutral oxygen efficiently cools the gas at temperatures in excess of a few 1000 K. With $n_{\rm O}$ denoting the neutral oxygen particle density the cooling rate is (Sternberg & Dalgarno 1989)

\begin{displaymath}\Lambda_{\rm OI~630nm} = 1.8 \times 10^{-24} n_{\rm O}~n_{\rm e}
\exp\left( -22~800/T_{\hspace*{-0.2ex}\rm g}\right) .
\end{displaymath} (107)

6.4 Miscellaneous heating/cooling processes

We list below two additional processes that can cause either heating or cooling of the gas.

6.4.1 Thermal accommodation

Following Burke & Hollenbach (1983), the energy exchange rate by inelastic collisions between grains and gas particles is

\begin{displaymath}\Gamma_{\rm g-g} = 4 \times 10^{-12}~\pi\langle a^2\rangle~n_...
...{\rm g}}~(T_{\hspace*{-0.2ex}\rm d}-T_{\hspace*{-0.2ex}\rm g})
\end{displaymath} (108)

For gas temperatures $T_{\hspace*{-0.2ex}\rm g}$ higher than dust temperatures $T_{\hspace*{-0.2ex}\rm d}$, this rate turns into a cooling rate $\Lambda_{\rm g-g}$. The thermal accommodation coefficient $\alpha_{\rm T}$ is set to the typical value for silicate and graphite dust of 0.3 (Burke & Hollenbach 1983).

6.4.2 free-free heating/cooling

Free-free transitions directly convert photon energy into thermal energy (ff-heating) or vice versa (ff-cooling) during electron encounters. The heating rate $\Gamma_{\rm ff}$ and cooling rate $\Lambda_{\rm ff}$ are given by
                                 $\displaystyle \Gamma_{\rm ff}$ = $\displaystyle 4\pi \int \kappa^{\rm ff}_\nu J_\nu~{\rm d}\nu$ (109)
$\displaystyle \Lambda_{\rm ff}$ = $\displaystyle 4\pi \int \kappa^{\rm ff}_\nu B_\nu(T_{\hspace*{-0.2ex}\rm g})\;{\rm d}\nu$ (110)
$\displaystyle \kappa^{\rm ff}_\nu$ = $\displaystyle n_{\rm e}^2~\sigma^{\rm ff}_\nu
+ n_{\rm e}~n_{\rm H}~\sigma^{\rm...
... H_2}~\sigma^{\rm H_2^-ff}_\nu
+ n_{\rm e}~n_{\rm He}~\sigma^{\rm He^-ff}_\nu ,$ (111)

where $\kappa^{\rm ff}_\nu$ is the free-free gas opacity [cm-1]. The free-free cross-sections $\sigma^{\rm ff}_\nu$ [cm5] for bremsstrahlung of singly ionised gases are taken from (Hummer 1988), for H-ff from (Stilley & Callaway 1970), for H2-ff from (Somerville 1964), and for He-ff from (John 1994).

7 Sound speeds

After the chemistry (see Sect. 5) and the thermal gas energy balance (see Sect. 6) have been solved throughout the disk volume, all particle densities ni and the kinetic temperature of the gas  $T_{\hspace*{-0.2ex}\rm g}$ are known, and P ROD IM O can update the isothermal sound speeds on the numerical grid cT2(rj,zk)as preparation for the next iteration of the hydrostatic disk structure (see Sect. 3).

                     $\displaystyle \rho$ = $\displaystyle n_{\rm e}~m_{\rm e} + \sum_i n_i~m_i$ (112)
p = $\displaystyle \Big( n_{\rm e} + \sum_i n_i\Big) ~kT_{\hspace*{-0.2ex}\rm g}$ (113)
cT2 = $\displaystyle p/\rho.$ (114)

gas in thermal balance & $T_{\hspace*{-0...
...ludegraphics[height=7.5cm,width=8.8cm]{1821fg10.eps} \end{tabular}
\end{figure} Figure 7:

Gas temperature structure $T_{\hspace *{-0.2ex}\rm g}(r,z)$ in a model for a T Tauri type disk with $M_{\rm disk}\!=\!0.01~M_\odot$ (l.h.s.). See further model parameter in Table 5. On the r.h.s. we show the results for the same parameter, if $T_{\hspace*{-0.2ex}\rm g}\!=\!T_{\hspace*{-0.2ex}\rm d}$ is assumed. The white dashed lines show $A_V\!=\!1$ and $A_V\!=\!10$.

Open with DEXTER

gas in thermal balance & $T_{\hspace*{-0...
...height=7.5cm,width=8.8cm]{1821fg12.eps} \end{tabular}\\ *[-3mm]
\end{figure} Figure 8:

Dust temperature structure $T_{\hspace *{-0.2ex}\rm d}(r,z)$. The differences between the l.h.s. and the r.h.s. model are caused by the different density structures which are a consistent result of the entire coupled physical problem.

Open with DEXTER

8 Results

We apply our P ROD IM O model to a typical passive protoplanetary disk of mass $M_{\rm disk}\!=\!0.01~M_\odot$ which extends from 0.5 AU to 500 AU. The central star is assumed to be a T Tauri-type ``young sun'' with parameters $T_{\rm eff}\!=\!5770~$K and $L_\star\!=\!1~L_\odot$, and to emit excess UV of predominantly chromospheric origin as shown in Fig. 2. The stellar UV excess creates an unshielded UV radiation strength of about $\chi\!=\!2\times10^6$ at 1 AU (see Eq. (41)). Further parameter of our model are summarised in Table 5. Our selection of elements and chemical species is outlined in Table 1, and the applied element abundances are listed in Table 3.

The model uses a $150 \times 150$ grid of points which are arranged along radial and vertical rays which enables us to calculate the respective column densities and line optical depths in a simple way. The spatial resolution is much higher in the inner regions and the grid points are also somewhat concentrated toward the midplane. About half of the grid points are located inside of 2.25 AU in this model to resolve the strong gradients in the radiation field and in the thermal and chemical structure occurring just inside of the inner rim.

Table 5:   Parameter of the model depicted in Figs. 7 to 14.

gas in thermal balance & $T_{\hspace*{-0...
...[height=7.5cm,width=8.8cm]{1821fg14.eps} \end{tabular}\\ *[-3mm]
\end{figure} Figure 9:

Density structure $n_{{\rm\langle H\rangle}}(r,z)$ as function of relative height z/r and $\log r$. The red dashed line on the l.h.s. encircles hot regions $T_{\hspace*{-0.2ex}\rm g}\!>\!1000~$K.

Open with DEXTER

8.1 Disk structure

The physical structure of the disk is a consistent result of all model components: dust radiative transfer, chemistry, and heating and cooling balance. In order to explore how important the inclusion of the gas heating and cooling balance is for the resulting disk structure, we compare the full model (depicted on the l.h.s. of the following figures) to a comparison model (r.h.s.) where we have assumed $T_{\hspace*{-0.2ex}\rm g}\!=\!T_{\hspace*{-0.2ex}\rm d}$ throughout the disk.

8.1.1 Thermal structure

Figures 7 and 8 show the resulting gas and dust temperature structures of the models, respectively. The most obvious feature in Fig. 7 is a hot surface layer ( $T_{\hspace*{-0.2ex}\rm g}\!\approx\!4000$-7000 K) which bends around the inner rim and continues radially to about 10 AU. This hot surface layer is situated above $z/r\!\ga\!0.13$ in this model. Its lower edge is not related to the vertical AV but rather to the position of the shadow casted by the puffed-up inner rim. It coincides with the first occurrence of CO and other molecules like OH (see Fig. 12). The hot surface layer is optically thin, predominantly atomic (molecule-free) and directly heated by the stellar radiation in various ways (see Sect. 8.1.4).

The shielded and cold regions in the midplane ( $z/r\!\la\!0.06$) are characterised by small deviations between  $T_{\hspace*{-0.2ex}\rm g}$ and $T_{\hspace*{-0.2ex}\rm d}$, due to effective thermal accommodation between gas and dust. However, beyond some critical radius, here $\approx$100 AU, even the midplane regions become optically thin, and the interstellar UV irradiation causes an increase of  $T_{\hspace*{-0.2ex}\rm g}$. We find midplane temperatures up to $T_{\hspace*{-0.2ex}\rm g}\!\approx\!2T_{\hspace*{-0.2ex}\rm d}$ around 400 AU in this model. The critical radius is related to $A_V\!=\!1$ and increases with disk mass.

The upper layers $z/r\!\ga\!0.1$ at $r\!\ga\!20~$AU show no clear trend, both $T_{\hspace*{-0.2ex}\rm g}\!<\!T_{\hspace*{-0.2ex}\rm d}$ and $T_{\hspace*{-0.2ex}\rm g}\!>\!T_{\hspace*{-0.2ex}\rm d}$ is possible, due to a complicated superposition of various heating and cooling processes.

Apart from the thermally decoupled layers at the inner rim, the surface and the very extended layers, the disk temperature is mainly controlled by the dust continuum radiative transfer (see Fig. 8). $T_{\hspace *{-0.2ex}\rm d}(r,z)$ shows all the features typical for protoplanetary disks (see e.g. Pinte et al. 2009; Pascucci et al. 2004). The midplane dust optical depth at $1~\mu$m is about $1.8\times 10^5$ in this model. The slightly different $T_{\hspace*{-0.2ex}\rm d}$-results for the two models are caused by the different density structures (see Fig. 9) which depend on  $T_{\hspace*{-0.2ex}\rm g}$. In case of the full model, the vertically extended inner regions scatter the star light and thereby heat the disk from above.

gas in thermal balance & $T_{\hspace*{-...
{1821fg16.eps} \end{tabular}
\end{figure} Figure 10:

Density structure of the puffed-up inner rim. Regions with $n_{{\rm\langle H\rangle}}\!<\!10^{~5.5}$ are suppressed for output (shown in white colour). The red dashed contour line on the l.h.s. encircles hot regions $T_{\hspace*{-0.2ex}\rm g}\!>\!1000~$K.

Open with DEXTER

8.1.2 To flare or not to flare

Figure 9 shows the resulting density structures of both models. The full model (l.h.s.) exhibits a remarkable vertical extension (up to $z/r\!\approx\!1$) of both the inner rim and the surface layers inward of $r\!\la\!10~$AU. According to Eq. (5), the vertical scale height H is approximately (assuming $z\!\ll\!r$, $c_T\!=~$const.) given by

\begin{displaymath}\left(\frac{H}{r}\right)^2 = \frac{2r~c_T^2}{GM_\star} ,
\end{displaymath} (115)

where H is defined as $\rho(z)\!\approx\!\rho(0)\exp(-z^2/H^2)$. The temperature ratio $T_{\hspace*{-0.2ex}\rm g}/T_{\hspace*{-0.2ex}\rm d}$ reaches values of about 10 - 30 in the hot inner rim and the surface regions, and since $c_T^2\!\propto\!T_{\hspace*{-0.2ex}\rm g}$, the disk is vertically more extended by about the same factor in comparison to the $T_{\hspace*{-0.2ex}\rm g}\!=\!T_{\hspace*{-0.2ex}\rm d}$-model. This applies to the inner rim in particular, because it is hot even at $z\!=\!0$, whereas the enhancing effect only starts at $z/r\!\ga\!0.1$ in general. However, in the regions $r\!\approx\!1$-7 AU, $T_{\hspace*{-0.2ex}\rm g}$is almost constant in the hot surface layer ($\approx$4000-7000 K) and so H/r increases further with increasing radius, and the disk reaches its maximum vertical extension here.

It is noteworthy that the vertical density structure $\rho(z)$ may be locally inverted. Since Eq. (5) is a pressure constraint, the density must locally re-increase if $T_{\hspace*{-0.2ex}\rm g}$ drops quickly with increasing height. This happens in the uppermost layers, in particular around 10 AU at $z/r\!\approx\!0.8$, a region which causes the most numerical problems during the course of the global iterations.

At larger radii $\ga$30 AU, both models show a comparable vertical extension, characterised by a generally flaring structure. The ``flaring'' (increase of H/r with increasing r) is a natural consequence of the radial dust temperature profile varying roughly like $T_{\hspace*{-0.2ex}\rm d}(r)\!\propto\!r^{-p}$ with $p\!\approx\!0.25$ in the midplane and $p\!\approx\!0.35~$- 0.45 in the optically thin parts, so $H/r\propto r^{~p}$.

8.1.3 The puffed-up inner rim

Figure 10 shows a magnification of the density structure in the innermost regions. The figure demonstrates the large impact of the treatment of the gas temperature in the model on the resulting disk structure. There is a rapid decline of the density between $n_{{\rm\langle H\rangle}}\!=\!10^9\rm ~cm^{-3}$ and $10^8\rm ~cm^{-3}$, which is caused by the steep $T_{\hspace*{-0.2ex}\rm g}$-increase at given pressure at the top of the shadow at $z/r\!\approx\!0.13$ casted by the inner rim (l.h.s.). Therefore, such densities merely exist in the model close to the star, but the cool and dense midplane regions $({>}10^9\rm ~cm^{-3})$ are surrounded by an extended ``halo'' composed of thin hot atomic gas of almost constant density ( $n_{{\rm\langle H\rangle}}\!=\!10^8$ to $10^7\rm ~cm^{-3}$) which extends as high up as $z/r\!\approx\!0.5$. These results are astonishingly robust against variation of the disk mass $M_{\rm disk}$ between 10-4 and $10^{-1}~M_\odot$ - we always find the same kind of halo composed of the same kind of gas with the same densities. Only the midplane regions contain more or less cold matter, according to $M_{\rm disk}$.

The assumed position of the inner rim at 0.5 AU in our model implies maximum dust temperatures of about 500 K, which is well below the dust sublimation temperature, and the shape of the inner rim is controlled by the radial force equilibrium at the inner edge which implies a smooth density gradient, see Sect. 3.1. In contrast, Isella & Natta (2005) investigated the effect of pressure-dependent sublimation of refractory grains on the shape of the inner rim. In reality, different kinds of refractory grains will be present which have not only different and pressure-dependent sublimation temperatures, but the dust temperatures are strongly dependent on dust kind due to dust opacity effects (see Woitke 2006), which can be expected to result in a highly complex chemical structure of the inner rim.

In comparison, the $T_{\hspace*{-0.2ex}\rm g}\!=\!T_{\hspace*{-0.2ex}\rm d}$-model does not possess the hot surface layers and, consequently, shows a much flatter structure. The inner rim is much less puffed-up causing the shadow borderline to be situated deeper. The inner ``soft edge'' is likewise less extended, only from 0.5-0.61 AU in the $T_{\hspace*{-0.2ex}\rm g}\!=\!T_{\hspace*{-0.2ex}\rm d}$-model, whereas is extends from 0.5-0.8 AU in the full model, or about 40% of the inner radius.

heating & cooling\\ [-1mm]
...egraphics[width=10.4cm,height=11cm]{1821fg18.eps} \end{tabular}
\end{figure} Figure 11:

Leading heating process (l.h.s.) and leading cooling process (r.h.s.) of the model in gas thermal balance. The black dashed contour line indicates an optical extinction of $A_V\!=\!10$.

Open with DEXTER

8.1.4 Thermal balance

Figure 11 shows the most important heating processes (l.h.s.) and the most important cooling process (r.h.s.) in the full model of the disk with the gas being in thermal balance. Again, there is a clear dividing line at $z/r\!=\!0.13$ coinciding with the shadow of the inner rim, which separates the directly illuminated hot surface layers from the shielded and cold midplane regions.

The central midplane of the disk below $A_V\!\approx\!10$ is dominated by thermal accommodation which assures $T_{\hspace*{-0.2ex}\rm g}\!\approx\!T_{\hspace*{-0.2ex}\rm d}$(see also Gorti & Hollenbach 2008; Kamp & Dullemond 2004; Nomura & Millar 2005). Since UV photons cannot penetrate these layers, cosmic-ray ionisation is the only remaining heating process, mostly compensated for by thermal accommodation cooling. In the central midplane $r\!\la\!1~$AU, before H2O freezes out (see Fig. 12), there is additionally H2O rotational cooling, as well as some H2quadrupole and CO rotational cooling just below $A_V \approx 10$.

Between $A_V\!\approx\!10$ and $z/r\!\approx\!0.13$, the UV radiation can partly penetrate the disk via scattering from above (see Fig. 4). This creates an active photon-dominated region with a rich molecular chemistry, where most of the abundant molecules like H2, CO, HCN, OH and H2O form, usually referred to as the ``intermediate warm molecular layer'' (Bergin et al. 2007). The layer is predominantly heated by H2 formation on grain surfaces and, with increasing height, by photo-effect on PAH molecules. The gas temperature increases upward in this layer, e.g. from $\sim$200 K to $\sim$700 K at 1 AU, but the additional heating can still be balanced by thermal accommodation in our model.

The upper edge of the warm molecular layer is characterised by a thin zone of intensive CO ro-vibrational cooling. Above this zone, CO is photo-dissociated - below this zone, the CO lines become optically thick. It is this CO ro-vibrational cooling that can counterbalance the upwards increasing UV heating for a while, until the heating becomes too strong even for CO. This happens just at the upper end of the disk shadow $z/r \approx 0.13$ where the direct stellar irradiation becomes dominant.

Above the CO layer, the temperature suddenly jumps to about 5000 K, all molecules are destroyed (thermally and radiatively), and we enter the hot surface layer described in the previous sections. This layer is predominantly heated by collisional de-excitation of vibrationally excited H$_2^\star $ (inner regions) and by PAH heating (outer regions). Although H2 is barely existent at these heights above the disk (concentration is 10-4 to 10-7, see Fig. 12), the few H2 molecules formed on grain surfaces can easily be excited by UV fluorescence, and these H$_2^\star $ particles undergo de-exciting collisions. This heating is balanced by various line cooling mechanisms. Since molecules are not available, atoms and ions like O  I and Fe  II are most effective. The non-LTE cooling by the wealth of fine-structure, semi-forbidden and permitted Fe  I and Fe  II lines has been investigated in detail by Woitke & Sedlmayr (1999), who found that in particular the semi-forbidden iron lines provide one of the most efficient cooling mechanisms for warm, predominantly atomic gases at densities $n_{{\rm\langle H\rangle}}\!=\!10^6$ to $10^{14}\rm ~cm^{-3}$.

Since the stellar optical to IR radiation can excite most of the Fe  II levels directly, radiative heating occurs. This ``background heating by Fe  II'' as referred to in Fig. 11 (l.h.s.) turns out to contribute significantly to the heating of the hot atomic layer close to the star ($r \la 10~$AU). In fact, further analysis shows that the gas temperature in a large fraction of the hot atomic layer is regulated by $\Gamma_{\rm FeII} \approx \Lambda_{\rm FeII}$, i.e. by radiative equilibrium of the gas with respect to the Fe  II line opacity. Similarly, we find a small zone in the midplane just behind the inner rim where radiative equilibrium with respect to the water line opacity is established. The regulation of the gas temperature via radiative equilibrium is a typical feature for dense gases in strong radiation fields, e.g. in stellar atmospheres. This behaviour is rather unusual in PDR and interstellar cloud research from where most of the other heating and cooling processes have been adopted.

The more distant regions $\ga$50 AU are characterised by an equilibrium between interstellar UV heating (photo-effect on PAHs) and [C II] $158~\mu$m, [O I] 63 and $144~\mu$m, CO rotational line cooling, and thermal accommodation (e.g.  Kamp & van Zadelhoff 2001).

8.2 Chemical structure

...e*{-6mm}\includegraphics[width=6cm]{1821fg30.eps} \end{tabular}
\end{figure} Figure 12:

Chemical composition of the gas in a T Tauri type protoplanetary disk with $M_{\rm disk}\!=\!0.01~M_\odot$, showing the concentrations $\epsilon_i\!=\!n_i/n_{\rm\langle H\rangle}$, where ni is the particle density of kind i and $n_{{\rm\langle H\rangle}}$ the total hydrogen nuclei density. Upper row: H, H2 and free electrons, second row: C+, C and CO, third row: O, OH and H2O and lower row: H2O ice, CO2 ice and CO ice. Note the different scaling for H, H2 and e- in the first row.

Open with DEXTER

The following discussion of the chemical results focuses on aspects that are relevant for an understanding of the two-dimensional disk structure. We restrict it to the most important atomic and molecular cooling species and the species that trace the dominant carriers of the abundant elements hydrogen, carbon and oxygen throughout the disk (Fig. 12). A more detailed discussion of particular chemical aspects and their relevance to observations will be the topic of future work.

8.2.1 Atomic and molecular hydrogen

Inward of about 10 AU, the H/H2 transition occurs at the lower boundary of the hot surface layer. There is a very sharp gradient of UV field and gas temperature explained by the shadow casted by the dust in the inner rim. Above the shadow, the gas temperature is high enough to efficiently destroy molecular hydrogen via $\rm H_2 + H \to 3H$, and also by collisions with atomic oxygen.

At larger distances, H2 can form on grain surfaces as soon as the dust temperature drops to about 100 K, where the formation efficiency $\epsilon(T_{\hspace*{-0.2ex}\rm d})$ increases sharply. This happens primarily in the secondary puffed-up regions around 10 AU. The formation of molecular hydrogen beyond this distance is mainly controlled by H2self-shielding, which is an intrinsically self-amplifying (i.e.  unstable) process. In addition, the gas density increases by a factor of $\sim$2 when H2 forms at given pressure, which causes increased collisional H2 formation rates in comparison to the photo-dissociation rates. This H2 formation instability leads to local overdense H2-rich regions in an otherwise atomic gas at high altitudes at about 10 AU in our model. Other molecules like OH and H2O are also affected and these molecules can show even larger concentration contrasts as compared to H2 which causes the instability.

8.2.2 Electron concentration and dead zone

The electron density in the upper part of the disk is set by the balance between UV ionisations and electron recombinations of atoms and molecules. In the UV obscured, cold and icy midplane below $z/r\!\approx\!0.05$, extending radially from just behind the inner rim to a distance of about 30 AU, the electron concentration drops to values below 10-8, but cosmic ray ionisations maintain a minimum electron concentration of $\sim$10-10 throughout the disk, because the vertical hydrogen column densities in this model are insufficient to absorb the cosmic rays ( $N_{\rm H_2}\!\approx\!10^{25}~\rm cm^{-2}$at $r\!=\!3~$AU). An electron concentration of $\sim$10-10 is two orders of magnitude larger than the minimum value of $\sim$10-12 required to sustain turbulence generation by magneto-rotational instability (MRI), see Sano & Stone (2002). Thus, our model does not possess a ``dead zone'' in the planet forming region, which is different from studies about massive and compact, actively accreting disks (e.g. Ilgner & Nelson 2006).

...ncludegraphics[width=9cm]{1821fg32.eps} \end{tabular}\\ *[-3mm]
\end{figure} Figure 13:

Cooling relaxation timescale $\tau _{\rm cool}$ (l.h.s.) and chemical relaxation timescale $\tau _{\rm chem}$ (r.h.s.).

Open with DEXTER

8.2.3 C$^+\!$, C, O, CO, OH, and H2O

Outside the shadowed regions, the models clearly show the classical C+/ C / CO / CO-ice transition as expected from PDR chemistry (e.g. Jonkheid et al. 2004; Gorti & Hollenbach 2008; Kamp & Dullemond 2004). However, there are some important differences to note in the 1-10 AU range. The dominant form of carbon in the midplane is CH4. At those high densities, oxygen is locked up into H2O-ice, leaving carbon to form methane instead of CO. Above the icy regions, in a belt up to $A_V \approx 10$, water molecules evaporate from the ice and CO becomes again the dominant carbon and oxygen carrier.

At radial distances between $\approx$0.8-10 AU in the warm intermediate layer, the model shows a double layer with high concentrations of neutral C, OH, H2O and other, partly organic molecules like CO2 HCN and H2CO (not depicted). This double layer is a result of the full 2D radiative transfer modelling in P ROD IM O . The radial UV intensities drop quickly by orders of magnitude at the position of the inner rim shadow ( $z/r \approx 0.13$). The UV radiation field then stays about constant, until $A_V\approx 1$ is reached, and also the vertical (+ scattered) UV intensities decrease. In combination with the downward decreasing gas temperatures and increasing gas densities, this produces two layers of hot and cold OH and H2O molecules with a maximum of C+ in between. van Zadelhoff et al. (2003) have undertaken similar investigations showing that dust scattering leads to a a deeper penetration and redirection of the stellar UV into the vertical direction, with strong impact on the photo-chemistry.

8.2.4 Ice formation

The ice formation is mainly a function of kind, gas density and dust temperature. Hence, the location of the individual ``ice lines'' strongly depend on the disk dust properties assumed, such as total grain surface area, disk shape and dust opacity. Water and CO2 ice formation is mostly restricted to the midplane, where the densities are in excess of 1010 cm-3, the reason being mainly the reaction pathways leading to the formation of the gaseous molecules that form these ice species. In addition, UV desorption counteracts the freeze-out of molecules in the upper layers at large distances from the star.

Inside 100 AU, densities are high enough to form water in the gas phase which subsequently freezes out onto the cold grains ( $T_{\hspace*{-0.2ex}\rm d}\la 100$ K). This is a consequence of our stationary chemistry that does not care about the intrinsically long timescale for ice formation (see Fig. 13). As densities drop and conditions for water formation in the gas phase become less favourable, oxygen predominantly forms CO, which freezes out at dust temperatures below $\sim$25 K at large distances. There is an intermediate density and temperature regime (20-100 AU), where significant amounts of CO2-ice are formed.

8.3 Timescales

An important question is whether our assumptions of gas energy balance and kinetic chemical equilibrium are valid in protoplanetary disks. The cooling relaxation timescale is calculated as

\begin{displaymath}\tau_{\rm cool} = \frac{3p}{2T_{\hspace*{-0.2ex}\rm g}}~\bigg...
...\partial Q}{\partial T_{\hspace*{-0.2ex}\rm g}}\bigg\vert^{-1}
\end{displaymath} (116)

where $Q\!=\!\Gamma\!-\!\Lambda$ is the net heating rate. The chemical relaxation timescale is more complicated. We calculate it as

\begin{displaymath}\tau_{\rm chem} = \max\limits_{{\rm valid}\ n}~\big\vert{\rm Re}\{\lambda_n\}^{-1}\big\vert
\end{displaymath} (117)

where $\lambda_n$ are the valid eigenvalues of the chemical Jacobian $\partial F_i/\partial n_j$ (see Eq. 64, without element conservation). Since F obeys $N_{\rm el}$ auxiliary conditions (element conservation), there are  $N_{\rm el}$ redundant modes which have (mathematically) zero eigenvalues. Numerics yields extremely small $\lambda$ for these modes, which must be disregarded.

Figure 13 (l.h.s.) shows that the cooling timescale in the disk is smaller than typical evolution timescales by orders of magnitude, and also smaller than typical mixing timescales (e.g. Ilgner et al. 2004), justifying our assumption of thermal balance. In particular, the gas is thermally tightly coupled to the dust via thermal accommodation in the midplane regions where the cooling timescale scales as $\tau_{\rm cool} \approx 1/\rho$.

The chemical relaxation timescales (r.h.s. of Fig. 13) show that apart from the icy midplane, where $\tau _{\rm chem}$ can be as large as 108 yrs, the chemical relaxation timescale is typically 104 yrs or shorter, and as short as $\approx$1-100 yrs in the photon-dominated warm intermediate layer, where most spectral lines form.

\end{tabular}\\ *[-3mm]
\end{figure} Figure 14:

Column density averaged gas temperatures $\langle T_{\hspace*{-0.2ex}\rm g}\rangle$ (Eq. (120), full lines) over vertical species column density for different radii labelled by the numbers in AU. The thin dashed lines show the column density averaged dust temperatures $\langle T_{\hspace*{-0.2ex}\rm d}\rangle$ for comparison. The small open circles, large open circles and large full circles mark vertical $A_V\!=\!0.1$, $A_V\!=\!1$ and $A_V\!=\!10$, respectively. The vertical arrows at the top of the figures mark the species column densities $N_{\rm thick}$ where the indicated lines become optically thick (Eq. (118)). We add some horizontal (rightward) arrows for rotational and ro-vibrational lines to demonstrate that $N_{\rm thick}$ is larger for high excitation lines. The leftward arrows on the r.h.s. of the plots indicate temperatures sufficient to thermally excite the upper level of the transitions $5~kT\!\approx\!E_u$.

Open with DEXTER

8.4 Gas emission lines

In order to discuss from which part of the disk the various gas emission lines come from, we provide some simple estimates of column densities and excitation temperatures in this section. Full 2D non-LTE line transfer calculations will be covered in subsequent papers.

The species column density $N_{\rm thick}$ required to achieve unit line optical depth can be calculated from Eq. (75). Assuming maximum and vanishing population in the lower and upper level, respectively ( $n_l\!\approx\!n_{\rm sp}$, $n_u\!\approx\!0$), the result for $\tau^{\rm ver}_{\rm ul}\!=\!1$ is

\begin{displaymath}N_{\rm thick} = \frac{g_l}{g_u}\frac{8\pi\nu^3\Delta v_D}{A_{\rm ul}~c^3}\cdot
\end{displaymath} (118)

Table 6 shows some typical values of $N_{\rm thick}$for permitted atomic resonance-lines, atomic fine-structure lines, and some rotational and ro-vibrational molecular lines. The actual critical column density is higher by a factor of $n_{\rm sp}/n_l$ if the population of the lower level is less than maximum.

Since the emission lines get saturated around $\tau\!\approx\!1$, the majority of the observable line flux originates in a surface region of thickness $N^{\rm ver}_{\rm sp}\!\approx\!N_{\rm thick}$. From the full model, we calculate the vertical species column densities $N^{\rm
ver}_{\rm sp}$ and column density averaged gas temperatures $\langle
T_{\hspace*{-0.2ex}\rm g}\rangle_{\rm sp}$ defined as

    $\displaystyle N^{\rm ver}_{\rm sp}(r,z) =
\int\limits_z^{z_{\rm max}(r)}\!\!n_{\rm sp}(r,z)~{\rm d}z$ (119)
    $\displaystyle \langle T_{\hspace*{-0.2ex}\rm g}\rangle_{\rm sp}(r,N^{\rm ver}_{...
...s_z^{z_{\rm max}(r)}\!\!T_{\hspace*{-0.2ex}\rm g}(r,z)~n_{\rm sp}(r,z)~{\rm d}z$ (120)

Similar to Eq. (120), we define the column density averaged dust temperature $\langle T_{\hspace*{-0.2ex}\rm d}\rangle_{\rm sp}$. The results are shown in Fig. 14. Treating the column above $N^{\rm
ver}_{\rm sp}\!=\!N_{\rm thick}$ as optically thin, ignoring deeper layers for the line formation, and assuming LTE, the column density averaged gas temperatures $\langle
T_{\hspace*{-0.2ex}\rm g}\rangle_{\rm sp}$ at depth $N^{\rm
ver}_{\rm sp}\!=\!N_{\rm thick}$ provide an estimate of the expected excitation temperature of the observable line flux.

Table 6:   Species masses in the disk and emission line characteristics. A velocity width of $\Delta v_D\!=\!1$ km s-1 is assumed.

CII: the most simple case in Fig. 14 is the 157.7$~\mu$m fine-structure line of C II. The C+ column density reaches a maximum value of $\sim$ $10^{17}\rm ~cm^2$, almost independent of radius r. This value falls just short of $N_{\rm thick}$, meaning that the line remains mostly optically thin throughout the entire disk. The column of emitting C+ gas is situated well above the optically thick dust in the midplane as indicated by the missing AV marks in Fig. 14. Since the outer regions have a much larger surface area as compared to the close regions, and the [C II] 157.7$~\mu$m fine-structure line can easily be excited even in a cold low density gas out to 500 AU ( $E_u~{\rm [K]}\!=\!91$, $n_{\rm
cr} \approx 3\times 10^3~$cm-3), the line flux is expected to be dominated by the outer regions.   The [C II] 157.7$~\mu$m line probes the upper flared surface layers of the outer disk.

OI: the column of atomic oxygen gas responsible for the [O I] 63.2$~\mu$m fine-structure line extends deeper into the midplane and reaches column densities of $\sim$ $10^{19}\rm ~cm^{-2}$, i.e. well beyond $N_{\rm thick}$. Therefore, the line can be expected to be optically thick in most cases. The line saturates before a dust visual extinction of $A_V\!=\!0.1$ is reached, i.e. also this line mainly probes the conditions above the optically thick midplane. Due to the higher energy Eu required to excite the upper level, the line is expected to originate mainly in regions inward of about 100 AU. We estimate this radius roughly by the requirement that columns must provide temperatures $5~kT_{\hspace*{-0.2ex}\rm g}\!\ga\!E_u$(see Fig. 14, arrows on r.h.s.). At a column density of $N^{\rm
ver}_{\rm sp}\!=\!N_{\rm thick}$, the line intensity adopts an excitation temperature of $\approx$40-70 K ($r \ga 30~$AU) and $\approx$ 500-1000 K ($r \la 10~$AU). In these inner regions, the [O I] 63.2$~\mu$m line partly forms in the hot atomic surface layer of the disk (see Sect. 8.1.1) where $\langle T_{\hspace*{-0.2ex}\rm g}\rangle_{\rm O}$ peaks to about 1000-4000 K at column densities $\sim$ $10^{16}...10^{18}\rm ~cm^{-2}$, although eventually the excitation temperature lowers again as the column of oxygen gas extends into deeper, cooler layers. The final value of $\approx$ 500-1000 K is $\approx$5-10 times higher than expected from $T_{\hspace*{-0.2ex}\rm g}= T_{\hspace*{-0.2ex}\rm d}$ models (see dashed lines in Fig. 14), stressing the necessity to include the gas energy balance in models for spectral interpretation. The [O I] 63.2$~\mu$m line forms in the thermally decoupled surface layers inward of about 100 AU, above $A_V \approx 0.1$.

CO: both the lowest rotational line $1\!\to\!0$ of CO as well as the CO fundamental ro-vibrational line $(v,J) = (1,1) \to (0,0)$ need a column density of about 10 $^{15}\rm ~cm^{-2}$ to saturate. These CO column densities are exceeded by about 2.5-4.5 orders of magnitude at all radii, i.e.  these lines are optically thick throughout the entire disk. However, high-J rotational or high-J ro-vibrational transitions test deeper layers because $n_l\!\ll\!n_{\rm CO}$ and hence $N_{\rm
thick}\!\gg\!10^{15}\rm ~cm^{-2}$. The observable line intensities are triggered mainly by the excitation criterion. For ro-vibrational transitions, for example, temperatures of the order of a 1000 K are only available inward of about 10 AU, where $\langle T_{\hspace*{-0.2ex}\rm g}\rangle$ has a strong negative gradient, i.e. the emitting CO is situated just below the borderline between the hot atomic and the warm intermediate disk layer, and we find values of about $\langle T_{\hspace*{-0.2ex}\rm g}\rangle_{\rm CO}(N^{\rm
ver}_{\rm CO}\!=\!N_{\rm thick}) \approx 1000~$K. In the more distant columns, the ro-vibrational lines are also optically thick but the line source function $S^{\rm L}_{\rm ul}$ is very small. The lower rotational lines until $5\!\to\!4$ originate in the entire disk, but for higher rotational transitions, the emitting region shrinks remarkably. The CO rotational and ro-vibrational lines are mostly optically thick and probe very different regions in the disk, depending on the excitation energy Eu.

H2O: the situation for the rotational water lines is quite different. Large amounts of H2O only form in deep layers, typically below $A_V\!=\!1$. However, between 1 AU and 100 AU, the total H2O column densities are limited by about $10^{14}{-}10^{16}\rm ~cm^{-2}$ because of ice formation in even deeper layers. Since the rotational water lines have larger $A_{\rm ul}$ in general and get optically thick sooner, these column densities seem sufficient to saturate the low-lying rotational lines out to about 30 AU. Because gas and dust temperatures are coupled in the deep layers, the excitation temperatures of the rotational water lines in LTE with the local dust temperatures. Therefore, the line-to-continuum ratio might be a problem for observations.

Higher rotational lines are difficult to excite and the origin of such lines is probably limited to regions $r \la 1~$AU. In particular, the inside of the inner rim is full of water and the temperatures here are high enough to excite a wealth of high-excitation rotational and probably also ro-vibrational lines, $\langle T_{\hspace*{-0.2ex}\rm g}\rangle_{\rm
H_2O}(N^{\rm ver}_{\rm H_2O}\!=\!N_{\rm thick}) \approx 200~{\rm
K} {-} 2000~$K between $r\!=\!0.6~$AU and r = 0.8 AU. Another interesting region is the vertically extended zone around 10 AU (see Figs. 9 and 12) which contains some amounts of hot water. We are currently investigating the impact of this hot water layer on the rotational lines as observable with Herschel (Woitke et al. 2009). The rotational H2O lines probe the conditions in the midplane regions $A_V \approx 1{-}30$, the inside of the inner rim, and possibly a vertically extended region with hot water around $10~{\rm AU}$ according to this model.

9 Conclusions and outlook

This paper introduces a new code, P ROD IM O, to model the physical, chemical and thermal structure of protoplanetary disks. The strength of the new code lies in a fully coupled treatment of 2D dust continuum radiative transfer, gas phase and photo-chemistry, ice formation, heating ${\rm\hspace*{0.7ex}\&\hspace*{0.7ex}}$cooling balance, and the hydrostatic disk structure. In particular, we use the calculated radiation field as input for the photo-chemistry and as background continuum for the non-LTE modelling of atoms, ions and molecules. The resulting gas temperatures determine the vertical disk extension, which in turn serves as input for the radiative transfer. Another advantage of the code is the robustness of its kinetic chemistry module which is applicable to densities between 102 to $10^{16}\rm ~cm^{-3}$; this makes possible to model complete disks ranging from about $\sim$0.5 AU to 500 AU.

Heating and cooling: the heating ${\rm\hspace*{0.7ex}\&\hspace*{0.7ex}}$cooling balance of the gas is the key to understand the vertical disk extension. The stellar UV irradiation, both direct and indirect via scattering, produces hot surface layers, in particular inside of $\approx$10 AU, even without X-rays. We have included large non-LTE systems for Fe  II and CO ro-vibrational line transitions to counterbalance this heating in the disk regions with high $\rho$ and high $T_{\hspace*{-0.2ex}\rm g}$, but find that these transitions also open new channels of radiative heating by the stellar optical - IR irradiation. The disk surface layers close to the star are in fact often stellar-atmosphere-like and characterised by radiative equilibrium. More work is required to identify the important heating ${\rm\hspace*{0.7ex}\&\hspace*{0.7ex}}$cooling processes in these layers and the gas inside of the inner dust rim.

Puffed-up inner rim and atomic halo: applying the new concept of ``soft edges'' to the inner rim of a T Tauri disk, the models show a highly puffed-up inner rim extending up to $z/r \approx 0.7$, and an extended layer of hot ( $T_{\hspace*{-0.2ex}\rm g}\approx 5000~$K) and thin ( $n_{{\rm\langle H\rangle}}\approx 10^{~7}{-}10^{~8}\rm ~cm^{-3}$) atomic gas reaching up to $z/r \approx 0.5$ in the innermost 10 AU. This ``halo'' is located above the intermediate warm molecular layer that surrounds the compact, dense, cold and icy midplane. It seems questionable that the highly puffed-up structures are hydrodynamical stable and further investigations must show how these features are related to mixing, winds and gas removal (Alexander 2008; Woitke et al. 2008).

Scattering and photo-chemistry: the dust grains in the puffed-up inner rim and the halo scatter the stellar UV light back onto the disk surface, which enhances the photo-chemistry and the photo-desorption of ice at larger distances.

Chemistry: the surface regions of the model reveal the classical PDR structure for H2, H, C+, C and CO. However, due to the full 2D UV radiation transfer in P ROD IM O, a complicated multi-layered structure results for H2O and other organic molecules like CO2 and HCN, which depend sensitively on the model parameters. The UV radiation field in the intermediate warm layer is reduced in two steps, first the puffed-up inner rim blocks the direct path of the radial photons from the star, and second the vertical and scattered photons are absorbed in the deeper layers. We find in particular two layers of hot and cold water molecules. Mixing by hydrodynamical motions is likely to smooth out such structures (Tscharnuter & Gail 2007; Ilgner et al. 2004; Semenov et al. 2006).

Cooling and chemical timescales: from the calculated relaxation timescales we conclude that the assumption of gas thermal balance and kinetic chemical equilibrium should be sufficient for the interpretation of most gas emission lines, although mixing may play a role. In the midplane, the chemical timescale is as long as 108 yr due to ice formation, but the spectral lines form predominantly in higher layers where $\tau_{\rm
cool}\!\approx\!1~...~100~$yr and $\tau_{\rm
chem}\!\approx\!1~...~10^4~$yr, depending on distance.

Emission lines: disk emission lines originate mostly in the thermally decoupled surface layers, where $T_{\hspace*{-0.2ex}\rm g}> T_{\hspace*{-0.2ex}\rm d}$. A simple analysis of the line characteristics and column densities shows that the [C II] 157.7$~\mu$m fine-structure line probes the outer flared surface of the disk, whereas the [O I] 63.2$~\mu$m line forms in slightly lower layers, but also probes the hot atomic gas inside of 10 AU. We identify this line as the key to test our models against observations. CO and H2O rotational lines probe the conditions in the intermediate warm molecular layer and, for low excitation lines, the outer regions. The origin of most H2O lines is possibly restricted to regions $\la$30 AU, partly because the water abundances are higher there, and partly because the lines are thermally excited only in these layers.

Observations: we intend to apply P ROD IM O to a large sample of observational disk emission line data that will be collected by the P ACS spectrometer on the HERSCHEL satellite, open time Key Program G ASPS. Based on calculated chemical and thermal disk structures, detailed non-LTE line radiative transfer calculations for the O  I and C  II fine-structure lines as well as some CO and H2O molecular lines will be carried out for analysis. We expect to be able to determine the gas mass in disks from the line data, and to find tracers for hot inner layers. In the next decade, P ROD IM O can be used to interpret A LMA data which probe the physical conditions and the chemical composition of the gas in the planet forming regions of protoplanetary disks.



... SUPA[*]
The Scottish Universities Physics Alliance.
... HOENIX-model[*]
... wrong[*]
For species that can be ionised with visual light like H-, it might be actually better to use our AV scale. But most photo-reactions occur in the UV and AV is just used as an auxiliary variable.
...(D'Alessio et al. 1998)[*]
Without further adjustments, the viscous heating rate according to Eq. (104) scales as $\Gamma_{\rm vis}\!\propto\!p$ at given radius r. Since all known cooling rates scale as $\Lambda\!\propto\!\rho^2$ in the low density limit, there is always a critical height z above which the viscous heating would dominate the energy balance and lead to ever increasing  $T_{\hspace*{-0.2ex}\rm g}$ (well above 20 000 K) with increasing height z. We consider this behaviour as an artefact of the concept of viscous heating and/or $\alpha $-viscosity.

All Tables

Table 1:   Elements and chemical species.

Table 2:   Adsorption energies and photo-desorption yields.

Table 3:   Assumed element abundances in (gas + ice)

Table 4:   Non-LTE model atoms, ions and molecules.

Table 5:   Parameter of the model depicted in Figs. 7 to 14.

Table 6:   Species masses in the disk and emission line characteristics. A velocity width of $\Delta v_D\!=\!1$ km s-1 is assumed.

All Figures

\end{figure} Figure 1:

Concept of global iterations in P ROD IM O. The circular arrows on the r.h.s. indicate sub-iterations. For example, the dust temperature structure needs to be iterated in the continuum radiative transfer.

Open with DEXTER
In the text

\end{figure} Figure 2:

Incident stellar intensity compiled from two sources: a P HOENIX solar model spectrum with $T_{\rm
eff}\!=\!5800~$K, $\log~g\!=\!4.5$, $Z\!=\!1$ (black line) and the chromospheric flux of ``young sun'' HD 129 333 (Dorren & Guinan 1994, red line). Note that the ionising and photo-dissociating flux is assumed to be restricted to the interval [91.2 nm,205 nm] which includes Ly$\alpha $ at 121.6 nm.

Open with DEXTER
In the text

\end{figure} Figure 3:

Benchmark for the dust continuum radiative transfer part. Vertical cuts of the calculated dust temperature structure $T_{\hspace *{-0.2ex}\rm d}(r,z)$ are shown at different radii of a disk with midplane optical depth $\tau_{\rm 1\mu{\rm m}}\!=\!10^5$ at $1~\mu$m. The P ROD IM O results are shown with blue diamonds in comparison to the MCFOST results (Pinte et al. 2006) with black lines.

Open with DEXTER
In the text

...ludegraphics[width=6.4cm]{1821fg06.eps} \end{tabular}\\ *[-3mm]
\end{figure} Figure 4:

Comparison of the UV radiation field strengths $\chi $ (Eq. (41)) between the simplified geometrical approach (l.h.s.), taking into account only vertical and radial dust extinction (see Eq. (47)), and the results of the full 2D radiative transfer including scattering (middle) for the $M_{\rm disk}\!=\!10^{-2}~M_\odot$ model. The two dashed contour lines show $A_V\!=\!1$ and $A_V\!=\!10$ as calculated from the minimum of the radial and vertical column densities. The r.h.s. figure shows the ratio  $\chi /\chi ^{\rm geo}$, indicating that the regions mostly affected by scattering are situated roughly between $A_V\!=\!1$ and $A_V\!=\!10$.

Open with DEXTER
In the text

\end{figure} Figure 5:

Different pumping and escape probabilities according to the predominantly radial irradiation and the predominantly vertical escape.

Open with DEXTER
In the text

\end{figure} Figure 6:

Continuum mean intensities as input for non-LTE modelling. The calculated band-mean mean intensities are shown for one particular point (r,z) 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. The red line shows the band-mean incident stellar intensities $\nu I_\nu ^\star $ and the blue line shows the incident interstellar intensities $\nu I_\nu ^{\rm ISM}$. The radiation field has two major components, the dust attenuated UV - near IR part, originating mainly from the star, and the self-generated mid - far IR part, originating from thermal dust emission in the disk.

Open with DEXTER
In the text

gas in thermal balance & $T_{\hspace*{-0...
...ludegraphics[height=7.5cm,width=8.8cm]{1821fg10.eps} \end{tabular}
\end{figure} Figure 7:

Gas temperature structure $T_{\hspace *{-0.2ex}\rm g}(r,z)$ in a model for a T Tauri type disk with $M_{\rm disk}\!=\!0.01~M_\odot$ (l.h.s.). See further model parameter in Table 5. On the r.h.s. we show the results for the same parameter, if $T_{\hspace*{-0.2ex}\rm g}\!=\!T_{\hspace*{-0.2ex}\rm d}$ is assumed. The white dashed lines show $A_V\!=\!1$ and $A_V\!=\!10$.

Open with DEXTER
In the text

gas in thermal balance & $T_{\hspace*{-0...
...height=7.5cm,width=8.8cm]{1821fg12.eps} \end{tabular}\\ *[-3mm]
\end{figure} Figure 8:

Dust temperature structure $T_{\hspace *{-0.2ex}\rm d}(r,z)$. The differences between the l.h.s. and the r.h.s. model are caused by the different density structures which are a consistent result of the entire coupled physical problem.

Open with DEXTER
In the text

gas in thermal balance & $T_{\hspace*{-0...
...[height=7.5cm,width=8.8cm]{1821fg14.eps} \end{tabular}\\ *[-3mm]
\end{figure} Figure 9:

Density structure $n_{{\rm\langle H\rangle}}(r,z)$ as function of relative height z/r and $\log r$. The red dashed line on the l.h.s. encircles hot regions $T_{\hspace*{-0.2ex}\rm g}\!>\!1000~$K.

Open with DEXTER
In the text

gas in thermal balance & $T_{\hspace*{-...
{1821fg16.eps} \end{tabular}
\end{figure} Figure 10:

Density structure of the puffed-up inner rim. Regions with $n_{{\rm\langle H\rangle}}\!<\!10^{~5.5}$ are suppressed for output (shown in white colour). The red dashed contour line on the l.h.s. encircles hot regions $T_{\hspace*{-0.2ex}\rm g}\!>\!1000~$K.

Open with DEXTER
In the text

heating & cooling\\ [-1mm]
...egraphics[width=10.4cm,height=11cm]{1821fg18.eps} \end{tabular}
\end{figure} Figure 11:

Leading heating process (l.h.s.) and leading cooling process (r.h.s.) of the model in gas thermal balance. The black dashed contour line indicates an optical extinction of $A_V\!=\!10$.

Open with DEXTER
In the text

...e*{-6mm}\includegraphics[width=6cm]{1821fg30.eps} \end{tabular}
\end{figure} Figure 12:

Chemical composition of the gas in a T Tauri type protoplanetary disk with $M_{\rm disk}\!=\!0.01~M_\odot$, showing the concentrations $\epsilon_i\!=\!n_i/n_{\rm\langle H\rangle}$, where ni is the particle density of kind i and $n_{{\rm\langle H\rangle}}$ the total hydrogen nuclei density. Upper row: H, H2 and free electrons, second row: C+, C and CO, third row: O, OH and H2O and lower row: H2O ice, CO2 ice and CO ice. Note the different scaling for H, H2 and e- in the first row.

Open with DEXTER
In the text

...ncludegraphics[width=9cm]{1821fg32.eps} \end{tabular}\\ *[-3mm]
\end{figure} Figure 13:

Cooling relaxation timescale $\tau _{\rm cool}$ (l.h.s.) and chemical relaxation timescale $\tau _{\rm chem}$ (r.h.s.).

Open with DEXTER
In the text

\end{tabular}\\ *[-3mm]
\end{figure} Figure 14:

Column density averaged gas temperatures $\langle T_{\hspace*{-0.2ex}\rm g}\rangle$ (Eq. (120), full lines) over vertical species column density for different radii labelled by the numbers in AU. The thin dashed lines show the column density averaged dust temperatures $\langle T_{\hspace*{-0.2ex}\rm d}\rangle$ for comparison. The small open circles, large open circles and large full circles mark vertical $A_V\!=\!0.1$, $A_V\!=\!1$ and $A_V\!=\!10$, respectively. The vertical arrows at the top of the figures mark the species column densities $N_{\rm thick}$ where the indicated lines become optically thick (Eq. (118)). We add some horizontal (rightward) arrows for rotational and ro-vibrational lines to demonstrate that $N_{\rm thick}$ is larger for high excitation lines. The leftward arrows on the r.h.s. of the plots indicate temperatures sufficient to thermally excite the upper level of the transitions $5~kT\!\approx\!E_u$.

Open with DEXTER
In the text

Copyright ESO 2009

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.