A&A 415, 643-659 (2004)
DOI: 10.1051/0004-6361:20034061

Transport processes and chemical evolution in steady accretion disk flows

M. Ilgner¹ - Th. Henning² - A. J. Markwick³ - T. J. Millar³

1 - Computational Physics, Auf der Morgenstelle 10, 72076 Tübingen, Germany
2 - Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
3 - Department of Physics, UMIST, PO Box 88, Manchester M60 1QD, UK

Received 9 July 2003 / Accepted 12 November 2003

Abstract
We study the influence of mass transport processes on the chemical evolution in a protoplanetary accretion disk. Local transport processes by advection as well as global transport processes by diffusion are taken into account. Concerning the multi-component system only diffusion in the vertical direction was taken into account. Depending on the transport properties, different schemes are explored to couple/decouple the physical and chemical evolution. Our model is based on a simplified description of hydrodynamics in terms of a steady 1+1-D-$\alpha$-disk model and includes the kinetics of an extended chemical network of about 250 species. We restrict our calculations to the inner planet formation zone within a distance to the central star of 10 AU.
Vertical mixing does change the global chemical evolution as it is demonstrated in detail through a discussion of the chemistry of sulphur-bearing molecules. In addition, the influence of the local disk structure on the chemical evolution is analysed. Our results demonstrate that the chemical evolution is influenced globally by mass transport processes. However, in addition to mass transport processes, information about the local conditions, which determine the kinetics, is still needed in order to analyze the chemical evolution.

Key words: solar system: formation - stars: circumstellar matter - ISM: molecules - ISM: abundances - ISM: evolution - protoplanetary disks

1 Introduction

It is difficult to give a quantitative comparision between the previous models describing chemical evolution in protoplanetary accretion disks; the chemical models are full of different intentions. For instance, the models of the Heidelberg group (Bauer et al. 1997; Finocchi & Gail 1997; Finocchi et al. 1997; Gail 1998, 2001) include treatments of ice mantle evaporation, dust evaporation, three-body reactions, ionisations by cosmic rays, radionuclides and UV photons. Later, in addition to annealing and carbon combustion the model of Gail (2002) accounts for carbon dust oxidation.

The approach used by Willacy et al. (1998) and Aikawa et al. (1999) is similar to Bauer et al. but does not account for dust destruction processes. Later, Aikawa & Herbst (1999, 2001) and Aikawa et al. (2002) added deuterium chemistry to their chemical model. The influence of the ionisation processes such as X-rays on the chemistry were introduced by Aikawa & Herbst (1999), and were studied in more detail by Markwick et al. (2002). Furthermore, Willacy & Langer (2000) discussed the importance of photoprocessing.

However, summarizing the improvements of recent research the models mainly focused attention on upgrading and extending the chemical networks as well as the implementation of additional kinetic processes.

A thorough treatment of physical processes within protoplanetary disks in parallel to chemical simulations was not the focus of many previous studies. This kind of imbalance was finally broken by Willacy et al. (1998) who introduced a 1+1-D-disk model in order to evaluate the chemical evolution based on that model. In a series of simulations, concerning the chemical evolution in accretion disks, the advective mass transport was taken into account, cf. Duschl et al. (1996), Bauer et al. (1997), Finocchi et al. (1997), Finocchi & Gail (1997), Willacy et al. (1998), Aikawa et al. (1999), Aikawa & Herbst (1999), and Markwick et al. (2002). However, vertical mass transport processes caused by diffusion were not considered.

This is the first of a series of papers aimed at dealing with mass transport processes in accretion disks in the framework of chemical evolution. The fluid in the model is continuously injected at the outer edge of the disk and the flow (and chemistry) is followed until the lower radial bound of our calculation, 1 AU. At smaller radii, as the material accretes on to the central protostar, temperatures are very large and we would need to take into account the physical and chemical destruction of the dust grains. It is this steady injection of material which allows us to adopt a Langrangian description so that our derived chemical distributions represent the steady state distributions of molecules throughout the disk. While the present paper is dedicated to steady accretion disk flows, the case of a non-steady accretion disk flow will be considered in a companion paper (Ilgner & Henning, in preparation).

We briefly summarise the basic concepts of accretion disk physics and astrochemistry as far as these concepts are needed for understanding the model developed here. The outline of the paper is as follows. In Sects. 2 and 3 we explain the physics and kinetics governing the chemical evolution of accretion disks. In Sect. 4 different techniques are introduced to model the advection transport and the transport by diffusion in a reacting flow system. Finally, in Sect. 5, we discuss our results and summarise our conclusions in Sect. 6.

2 Disk physics

2.1 Transport in disks

The dynamics of accretion disks are closely related to the transport of energy and angular momentum. Although there is still a fundamental problem on the origin of this transport process (hydrodynamical vs. magneto-hydrodynamical origin), one has to have a very efficient transport mechanism such as turbulent motion. The idea that disks might be turbulent to produce the anomalous viscosity seems to be feasible since Balbus & Hawley (1991) have given an explanation for the origin of turbulence via the magneto-rotational instability (MRI). Nevertheless, MRI cannot work throughout the entire since there are "dead zones" where other processes must operate.

If one wants to study the chemical evolution within the accretion flow one cannot incorporate the latest developments in disk modelling. Despite the improved computational resources, we are not able to perform 3D calculations taking into account a complete network of chemical reactions. Thus, we are forced to simplify the description of disk physics including the mass transport processes. For that reason, we have to consider the development of recent accretion disk theory in order to get an idea about how the physics might be simplified.

As soon as we move to a two-dimensional (2D) description we are not in a position to track comprehensively the causes which lead to turbulence because turbulence is a "rotational'' (vorticities), and a three-dimensional phenomenon. Of course, one can infer a high value of the phenomenological turbulent kinematic viscosity $\nu_{{\rm t}}$ due to the presence of a high Reynolds number flow which is typically for turbulent flows. In order to keep the transport properties of a turbulent flow in a 2D model one might be attempted to replace the kinematic viscosity $\nu$ with the turbulent kinematic viscosity $\nu_{{\rm t}}$ within the Navier Stokes equation. However, one has to remember that turbulence is a dynamical phenomenon and not a feature of the fluid itself. The viscous stress tensor $\sigma_{r\varphi}^{\prime}$ (cylindrical coordinates r, $\varphi$, and z are assumed) is approximated by means of the standard $\alpha$-prescription given by Shakura  & Sunyaev (1973), i.e., $\sigma_{r\varphi}^{\prime}\!=\!\alpha P$, where $\alpha$ and P denote a dimensionless quantity and the gas pressure, respectively. Finally, a turbulence model describing the turbulent kinematic viscosity $\nu_{{\rm t}}$ is explicitly imposed in the form $\nu_{{\rm t}} \! \sim \! \alpha c_{\rm s} H$ (Pringle 1981), where $c_{\rm s}$ and H denote the sound speed and a suitably defined disk height, respectively.

As stated above, two-dimensional models of accretion disks are qualitatively quite different from the three-dimensional models in terms of the physical interpretation of turbulence. Under the assumption of an axisymmetric (i.e., $\partial$/ $\partial \varphi =0$) hydrodynamical 2D description of an accretion disk, we are especially interested in the flow pattern. For instance the results of Kley et al. (1993), who studied the properties of convectively unstable disks, can be utilised to get general information about the direction of the flow. As long as Schwarzschild's criterion is satisfied, the matter inside a convective region sinks at low temperatures and rises at high temperatures, respectively. The values of the vertically integrated local mass flow are characterised by strong variations. In addition, the vertical stratification is nonsymmetric with respect to the disk midplane. The limiting case that the disk is kept stable against convection results in a laminar flow pattern. It has been demonstrated by Kley & Lin (1992) that the variation of the radial inflow velocity with height approaches a steady state. For such a system, the vertically integrated local mass flow reaches an equilibrium value and the radial velocities represent inflow.

Finally, one can describe the global transport processes within an accretion disk, taking into account only a few of the equations determining the 2D radiation-hydrodynamics. The geometry of the disk is considered in the "thin disk'' approximation, assuming that the thickness of the disk is everywhere small compared to its radius. Based on the equation of continuity, the z-component of the momentum equation, and the balance equation of angular momentum one can derive a parabolic differential equation determining the dynamical evolution of the disk (Lynden-Bell & Pringle 1974):

$$\frac{\partial \Sigma}{\partial t} = \frac{3}{r} \frac{\part... ...left( \nu_{{\rm t}}\Sigma r_{}^{\frac{1}{2}} \right) \right],$$ (1)

where $v_{\rm r}$ is given by

$$v_{\rm r} = -\frac{3}{\Sigma r_{}^{\frac{1}{2}}} \frac{\part... ...al r} \left( \nu_{{\rm t}} \Sigma r_{}^{\frac{1}{2}} \right).$$ (2)

The quantities $\Sigma$ and $v_{\rm r}$ denote the surface density and the radial drift velocity, respectively. The viscous shear stress is approximated by the $\alpha$-prescription of Shakura & Sunyaev (1973) and the motion in azimuthal direction is a Keplerian motion.

2.2 Steady 1+1-D-disk model

Provided that $\nu_{{\rm t}} \Sigma$ is a function of $\Sigma$ itself, the steady-state condition applied to Eq. (1) will only define the quantity $\nu_{{\rm t}} \Sigma$, not $\Sigma$. The steady-state behaviour of Eq. (1) is given by Lynden-Bell & Pringle (1974):

$$\nu_{{\rm t}} \Sigma := \frac{\dot{M}}{3\pi} \left[ 1- \left( \frac{R_{\star}}{r} \right)_{}^{\frac{1}{2}} \right],$$ (3)

where $R_{\star}$ refers to the radius of the star, assuming that the star rotates more slowly than Keplerian at the inner boundary of the disk. In the steady-state case the evolution is regulated by one parameter only, the mass accretion rate $\dot{M}$.

The mass accretion rate $\dot{M}$ is constant for each ring of a stationary, cylindrically symmetric Keplerian disk. Extraction of energy and angular momentum from test particles results in a spiraling inward motion and a constant mass transport towards the central object. The loss of angular momentum can be understood in terms of shearing stresses. The energy will dissipate, and the corresponding term in the dissipation function $\Phi$ is given by:

 

$$\Phi \nu_{{\rm t}} \Sigma r_{}^2 \mid \nabla \Omega \mid_{}^2$$ = (4)

$$\frac{3}{4\pi} \frac{G M_{\star} \dot{M} }{r_{}^3} \left[ 1 - \left( \frac{R_{\star}}{r} \right)_{}^{\frac{1}{2}} \right],$$ =

refering to dissipation of energy per unit time and per unit area. $\Omega$ denotes the Keplerian value of the angular velocity, G and $M_{\star}$ the gravitational constant and the mass of the central object, respectively. The dissipated energy will be radiated locally somewhere above the disk midplane. Assuming that the radiation field is close to the blackbody value, i.e. $\pi B(T)\!=\!\sigma_{{\rm R}}T_{{\rm eff}}^4$, one can calculate the radiation flux emergent from one side of the disk surface:

$$\sigma_{{\rm R}} T_{{\rm eff}}^4 = \frac{3}{8\pi} \frac{G M... ...- \left( \frac{R_{\star}}{r} \right)_{}^{\frac{1}{2}} \right],$$ (5)

where $\sigma_{{\rm R}}$ and B refer to the Stefan-Boltzmann constant and the integrated Planck function. In addition, the energy flux at the disk midplane ($z\!=\!0$) vanishes because of reflection symmetry.

The still unknown upper boundary $z_{}^{\ast}(r)$ can be determined by considering the vertical disk structure which is determined by a two-point boundary problem. In particular, the z-components of mass density $\varrho$, pressure P, temperature T, and the radiative flux of energy F of the steady solution of the equations of radiation hydrodynamics have to be calculated up to the disk surface:

   

$$\frac{{\rm d}P}{{\rm d}z} - \varrho g_z$$ = (6)

$$\frac{{\rm d}T}{{\rm d}z} \left\{ \begin{array}{l@{\quad {\rm for} \quad}l} -\frac{\overlin... ...\nabla_{{\rm conv}} & \nabla_{{\rm rad}} > \nabla_{{\rm ad}} \end{array}\right.$$ = (7)

$$\frac{{\rm d}F}{{\rm d}z} \frac{9}{4} \varrho \nu_{{\rm t}} \frac{G M_{\star}}{r_{}^3}\cdot$$ = (8)

The quantities g_z, $\lambda$, and $\overline{\kappa}_{{\rm R}}$ denote the z-component of the acceleration of gravity, the fluxlimiter, and the Rosseland mean opacity. Further we used a compilation given by Semenov et al. (2003), in particular composite dust aggregates with a "normal'' chemical composition refering to the iron abundance within the silicates, i.e., Fe/(Fe+Mg) $\!=\!0.3$, were considered. Energy transport by radiation as well as by convection is assumed depending on the local ratio between the radiative and adiabatic temperature gradients $\nabla_{{\rm rad}}$ and $\nabla_{{\rm ad}}$, respectively. Considering radiation transport, we used a flux limited diffusion approximation of the radiation transport given by Levermore & Pomraning (1981). Based on this theory the magnitude of the flux is regulated to be no greater than the density times the maximum transport speed. Finally, the system of equations can be closed with the equation of state of an ideal gas.

In order to solve the coupled system of three first-order ordinary differential Eqs. (6-8), three boundary conditions for $P(z_{}^{\ast})$, $T(z_{}^{\ast})$, and $F(z_{}^{\ast})$ are needed. While the value of F can be obtained from Eq. (5), the boundary values of P and T are linked to the value of the optical depth $\tau_{}^{\ast}$ at $z_{}^{\ast}$. In Appendix we summarize the radiative transfer problem and the $\tau_{}^{\ast} \!=\! 2/3$ condition. Concerning the boundary condition for T, we get $T(z_{}^{\ast})\!=\!T_{{\rm eff}}$. The boundary problem can be reduced to find the photospheric disk height $z_{}^{\ast}$ by shooting methods. More than 100 zones are needed to calculate the vertical disk structure.

The changes in radial structure of the disk are determined by the viscous time scale. This time scale is much longer than the time scales for other relevant physical processes in the disk. Because of the difference between the viscous time scale and the other time scales, one may decouple the radial and vertical direction of a 2D model resulting in a 1+1-D-description. We assume implicitly that a static disk atmosphere will be established instantaneously during the viscous evolution. Of course, the differences between a 2D and 1+1-D-description still remain. However, in the steady-state case, one can expect the 1+1-D-model to represent the 2D-model.

2.3 Model calculations

To keep the underlying assumptions of the disk model valid, one is forced to restrict the inner as well as the outer boundary of the computational domain $[r_{{\rm min}},r_{{\rm max}}]$. The disk has to be kept optically thick, i.e., $\tau(r,z) \ge 2/3$, and stable against self-gravitation. While the instability against self-gravitation can be measured by Toomre's parameter Q (Toomre 1964), the optical thickness at the upper boundary $z_{}^{\ast}$ is fixed by the boundary condition $\tau_{}^{\ast} \!=\! 2/3$. The regions where the disk becomes optically thin can be estimated by means of approximate relations shown in Bell et al. (1997). The quantity $r_{{\rm max}}$ was chosen to be 10 AU which does not necessarily mean an upper limit concerning $\tau$ and Q. Depending on the parameters of the system ($\alpha$, $\dot{M}$, opacity table) $r_{{\rm max}}$ can be extended by a factor of 2. We also note that the inner regions of the disks adjust to a quasi steady state more quickly than the outer parts. Thus, the steady-state assumption is easier fullfilled in the inner regions.

Purely for kinetic reasons the inner boundary $r_{{\rm min}}$ was fixed at 1 AU. Due to destruction of the dust particles by evaporation, the the opacity will change and the mathematical description of the kinetics, based on a system of ordinary differential equations (see Eq. (12)), is no longer suitable. It finally results in a system of differential algebraic equations. The only computations in this field have been carried out by the Heidelberg group, e.g., Bauer et al. (1997).

The focus of our disk model is not to improve existing disk models describing the vertical structure. In fact, we balanced the recent developments of more sophisticated disk models against simpler models in order to include the main features. Our disk model is quite similiar to that of Bell et al. (1997); differences are additional heating sources. Ionization by energetic particles and the radioactive decay of ²⁶Al were taken into account (see e.g. D'Alessio et al. 1998).

Figure 1: The radial distribution of the photospheric temperatures (dotted line) in comparison to the model by D'Alessio et al. (1998) - (dashed line - for $\dot{M}\!=\!10^{-8}\ M_{\odot}/{\rm yr}$, $\alpha\!=\!0.01$, and the opacity table derived by Pollack et al. (1994). In addition, the temperature along the "superheated'' surface layer and the disk interior temperature after Chiang & Goldreich (1997) are plotted.

Open with DEXTER

Our model differs from the model of D'Alessio et al. concerning the influence of the stellar radiation on the disk structure. In this study we did not consider stellar heating sources. When we compare the physical parameters of our non-irradiative disk model with that of the irradiative model of D'Alessio for the regions between 1 and 10 AU, we find that the differences in temperature, density and optical depth are not large (see Fig. 1) and are similiar to the differences which are caused by, for example, adopting a different grain model and therefore different opacity model. That is to say, in the inner disk interior to 10 AU, viscous dissipation rather than radiation dominates the heating. The largest differences can be found in the radial distribution of the photospheric temperature $T(z^{\ast})$, see Fig. 1. In addition, the photospheric height $z^{\ast}$ is always below the height z_s where the stellar radiation is deposited. The temperature inversion caused by the stellar radiation field occurs only above the photospheric height $z^{\ast}$. An inversion of temperature occurs in the simple disk model of Chiang & Goldreich (1997). Although this model does not account for transport processes, it completes the disk models used in astrochemical modelling of accretion disks, e.g., Aikawa et al. (2002) (D'Alessio disk model), Willacy & Langer (2000) (Chiang-Goldreich disk model). In the case of the radial temperature distribution shown in Fig. 1, we also plot the temperature distribution of the "superheated'' surface layer $T_{{\rm ds}}$ and the interior disk temperature $T_{{\rm int}}$ applying the description given by Chiang & Goldreich. The stellar parameters are the same as the corresponding parameters used in the dynamical disk models.

  
3 Thermodynamics

Chemical reactions as well as diffusion processes belong to the class of irreversible processes. They are macroscopically described by non-equilibrium thermodynamics. Irreversible transformations are distinguishable from reversible transformations by the state function S, the entropy of the system. A local description of S requires a continuous function of space coordinates, i.e., $S \! = \! \int_{}^{V} \varrho s {\rm d}V$, where $\varrho$ and s denote the mass density and the entropy per unit mass, respectively. In particular, the rate of entropy change inside the system, i.e., ${\rm d}_iS/{\rm d}t \! = \! \int_{}^{V} \sigma {\rm d}V$, where $\sigma$ denotes the entropy source strength or entropy production per unit volume and unit time, is related to irreversible processes inside the volume element V. The entropy production is always a non-negative quantity, i.e., $\sigma \ge 0$, and the entropy production vanishes only for reversible processes.

A macroscopic expression for $\sigma$ can be obtained by $\sigma~=~\sum_i J_iX_i$, where J_i and X_i are the Cartesian components of the independent thermodynamic fluxes and forces. The fluxes and forces are only up to 2nd tensorial order and, therefore, the entropy production $\sigma$ can be split into four different products of the fluxes and forces. The product of J_i and X_i can be determined by scalar quantities as well as by (polar and axial) vectors, and by two symmetric tensors with zero trace. However, the forces and fluxes of different tensorial character do not couple (Curie Principle). According to this classification, the entropy production $\sigma$ caused by chemical reactions is related to scalar phenomena in contrast to vectorial phenomena like diffusion.

In the following two subsections, we introduce the phenomenological equations in order to get information about the thermodynamic forces driving chemical reactions and diffusion, respectively. In a first approximation, the phenomenological equations can be written as linear relations between the thermodynamical forces and fluxes

$$J_i=\sum_{j=1}^{n} L_{ij} X_j,$$ (9)

where the components of the tensor L_ij are called phenomenological coefficients. Under the assumption of the existence of spatial symmetry properties one can demonstrate that the Cartesian components of the fluxes depend only on a few of the components of the forces (Curie principle).

3.1 Chemical reactions

Within thermodynamics, chemical reactions are attributed to a class of phenomenological equations which account for scalar fluxes and scalar forces. The scalar flux of a fixed chemical reaction is given by

$$J_i = -\frac{1}{T}\sum_{j=1}^{r}L_{ij}A_j -\frac{1}{3T}L_{iv} \nabla \cdot \vec {v} \quad (i=1,\dots,r),$$ (10)

where T is the temperature and A_j refers to the chemical affinity of the reaction j introduced by

$$A_i = \sum_{j=1}^{n}\mu_{j}\nu_{ji}.$$

The quantities $\mu_{j}$ denote the chemical potential of the component j; the quantity $\nu_{ji}$ divided by the molecular mass M_i is proportional to the stoichiometric coefficient of the component i in reaction j. By means of the general form of the phenomenological Eq. (9), one can specify the thermodynamical forces $X_i\!=\!-A_{i}/T$. The second term in Eq. (10) corresponds to cross-phenomena because the scalar flux may be caused by chemistry as well as by bulk viscosity. This phenomenon must be taken into account in the case of an incompressible fluid or a fluid with vanishing volume viscosity (Maxwell fluid).

The rate of change of the mass of component i is determined by the scalar thermodynamical flux. In terms of the molar concentration $N_i\!=\! \varrho_i/M_i$ (where $\varrho_i$ and M_i denote the mass density and the molecular mass of component i, respectively) the law of mass action is given by

$$\frac{\ \partial}{\partial t}N_i = \sum_{j=1}^r \overline{\nu}_{ij} {\overline{J}}_j \quad \quad (i=1,\dots,n),$$ (11)

neglecting any macroscopic transport. Because of the prefered description, using molar concentrations, the quantity $\overline{\nu}_{ij}$ is the stoichometric coefficient and ${\overline{J}}_j$ the rate of the chemical reaction j. With respect to chemical kinetics Eq. (11) can be written as

$$\frac{\ \partial}{\partial t} N_i = \sum_{j=1}^r \overline{\... ...k_j^{^-} \prod_{l=m+1}^n N_l^{ \overline{\nu}_{lj}} \right),$$ (12)

where k_j⁺ and k_j^- refer to the rate constants for the forward and backward reactions, respectively.

3.2 Diffusion

Diffusion is a vectorial phenomena. In general, the diffusion flow $\overline{\vec{J}}_i^a$ of the component i

$$\overline{\vec{J}}_i^a = N_i\left( \vec{v}_i-\vec{v}_{}^a \right)$$ (13)

is defined with respect to various reference velocities $\vec {v}_{}^a~=~\sum a_i\vec{v}_i$ where a_i refers to the normalized weights of the component velocities $\vec{v}_{i}$.

Only the most simple case will be considered here, neglecting external forces and vectorial cross-phenomena. The thermodynamical force is determined by $\vec{X}_i~=~\left. (\nabla \mu_i) \right\vert _{{\rm T,P}} / T$ where $\mu_i$ denotes the chemical potential of component i. Because of the linear relations between the fluxes and thermodynamical forces in Eq. (9), the diffusion flow of component i, which is given by

$$\overline{\vec{J}}_i^a = - \sum_{j=1}^{n-1} D_{ij} \nabla N_j,$$ (14)

does not only depend on the gradient of the chemical potential of component i. Therefore, we have to deal with a tensor D_ik of (n-1)² diffusion coefficients. Only in the special case of diffusion in a binary system ($n\!=\!2$), the tensor can be reduced to a scalar obtaining

$$\begin{eqnarray*}\overline{\vec{J}}_1^a & = & - D \ \nabla N_1 \\ D & \equiv & ... ...{a_2^2 \varrho_1 T} \frac{\partial \mu_1}{\partial N_1} \cdot \end{eqnarray*}$$

In addition, the relation between the diffusion coefficient D and the phenomenological coefficient L₁₁^a is shown in the last formula.

In order to achieve a simple model of diffusion in a multi-component system, we are interested in decoupling the thermodynamical forces and fluxes. Decoupling can be forced if the entire system is approximated by the superposition of n self-diffusion models. Every self-diffusion model combines one carrier material and only one tracer material assuming that both components can be characterized by the same dynamical properties. Finally, the change of mass of component i caused by diffusion results in

 

$$\frac{\partial \ }{\partial t} N_i - \nabla \cdot \overline{\vec{J}}_i^a$$ = (15)

$$\nabla \cdot \left( D \ \nabla N_i \right) ,$$ =

where D denotes the binary diffusion coefficient. However, this form of the diffusion law, conventionally used in the kinetic theory of gases, is slightly different from, e.g., the description which can be found in Xie et al. (1995). The different descriptions come from the chosen reference velocity $\vec{v}^a$ and the normalized weights a_i, respectively. The description given by Xie et al. refers to "relative diffusion flows'' where the reference velocity is measured with respect to the $n{{\rm th}}$ component.

Here we assume that the disk is in a turbulent state. We now have to specify the diffusion coefficient D. A viscous fluid will support stresses which are not only directed normal to the surface of the fluid element. Normal stresses as well as tangential stresses are represented by the viscous stress tensor ${\bf\sigma}_{ij}^{\prime}$. The viscous stress tensor is part of the momentum flux density tensor ${\rm\Pi}_{ij}$ which accounts for the rate of momentum flow through the boundary. Remembering the momentum equation, one may interpret viscous friction by an exchange of momentum. Because of the structure of the momentum equation one can draw an analogy between the exchange of momentum and a diffusion process where the kinematic viscosity $\nu$ acts as a diffusion coefficient.

The kinematic viscosity $\nu$ is now replaced by the turbulent kinematic viscosity $\nu_{{\rm t}}$. By means of the Schmidt number ${\rm Sc}~=~\nu_{{\rm t}} / D$, measuring quantitative difference between the transport of material and momentum, one can determine the diffusion coefficient of Eq. (15). For example, the Schmidt number is less than unity if the spreading rate of material is greater than that of momentum.

However, the diffusion coefficient D is influenced by dust particles mixed within the gas. Therefore, the coupling between the gas and dust particles implies that one should define an effective diffusion coefficient. However, recent studies by Schräpler & Henning (in preparation) have shown that this effect on the diffusivity becomes important only for dust sizes above $0.01 \ {\rm cm}$, depending on the local position within the disk ( $r \! \le \! 10 {\rm AU}$). Therefore, we can neglect the role of dust and take the Schmidt number in our calculations to be ${\rm Sc}\!=\!1$. Experimental data typically account for Schmidt numbers ${\rm Sc}\! \approx \! 0.7$ (cf. Spalding 1971); our estimate of ${\rm Sc}\!=\!1$ will underestimate the efficiency of turbulent material transport due to mixing. Previous calculations of mixing in disks have also adopted ${\rm Sc}\!=\!1$ (Gail 2001; Wehrstedt & Gail 2002).

The eddy diffusion profile can be obtained by the disk model and the turbulence model, respectively. Assuming the $\alpha$-turbulence model, the local value of the turbulent viscosity $\nu_{{\rm t}}\!=\!\alpha c_{\rm s} \min[H_{\rm P},H]$ is determined by the model itself where $c_{\rm s}$ denotes the sound speed, $H_{\rm P}$ the pressure scale height, and H the disk height. Therefore, a priori estimates of the turbulent viscosity are not needed, in contrast to the eddy diffusion profiles in molecular clouds, cf. Rawlings & Hartquist (1997).

4 Modelling of a reacting flow system

In order to evaluate the importance of transport processes on the chemical disk evolution, it may be useful to start with a few preliminary remarks. It is obvious that the modeling of chemical disk evolution within a protoplanetary disk becomes simpler if one excludes mass transport processes. Remarkably, since the first studies concerning the chemical disk evolution, this simplified description is still quite popular, although models including mass transport processes are already available, e.g., Duschl et al. (1996), Willacy et al. (1998), Aikawa et al. (1999). In addition, we have to point out that models without transport processes discussed by Aikawa & Herbst (1999, 2001) and Willacy & Langer (2000) refer to a very specific application of disk models. While Willacy and Langer applied the Chiang-Goldreich model, Aikawa and Herbst utilized the so-called minimum-mass solar nebula model (Hayashi et al. 1985). However, a protoplanetary disk has to include the accretion process onto the central object and, therefore, the accretion flow within the disk has be taken into account.

Transport processes can be distinguished into transport due to advection and transport due to diffusion. This kind of distinction can be introduced by applying Eq. (13) to the equation of continuity, which for component i is:

$$\frac{\partial \ }{\partial t} N_i + \nabla \cdot \left( N_i... ...}^a \right) = - \nabla \cdot \overline{\vec{J}}_i^a + R_i \ ,$$ (16)

where R_i denotes the rate term for production and destruction of component i, respectively. The source term on the left hand side, $\nabla \cdot \left( N_i \vec{v}_{}^a \right)$, refers to advective transport while diffusive transport is represented by the source term on the right hand side, $\nabla \cdot \overline{\vec{J}}_i^a$.

As far as the kinetics is determined by non-equilibrium thermodynamics, simulations of the chemical disk evolution have not yet been performed in the Eulerian description. Simulations, assuming chemical equilibrium, have been carried out in order to study the influence of annealing and carbon combustion on the disk evolution, see Gail (2001). In a recent study, Gail (2002) extended his model, to include a restricted chemical network of 18 species in the gas phase.
In this study only vertical mixing of chemical species is account for diffusion processes; diffusion in radial direction would need a different formulation than we used, and would require a fully two-dimensional reaction-diffusion model. However, such work is in progress.

4.1 Advection

In this subsection we focus only on the advective transport of chemical species radially. As long as the turbulent mixing of species is excluded, the computational domain can be splitted into a set of horizontal layers. Once vertical mixing is taken into account, the vertical grid is needed.

Instead of the Eulerian description, the Lagrangian description will be used here to study the time variation of the Lagrangian fluid element. The Eulerian description can be transformed into the Lagrangian description by means of the substantial derivative ( ${\rm d}/ {\rm d}t$) which is simply related to the time derivative ($\partial$/ $\partial t$) computed at fixed gridpoints in space via

$$\noindent \frac{{\rm d}\ }{{\rm d}t} = \frac{\partial \ }{\partial t} + \left( \vec{v} \nabla \right).$$

Thus, the advective term on the left hand side of Eq. (16) can be replaced by $N_i \nabla \cdot \vec{v}_{}^a$.

As yet, the coupling between a particular fluid parcel following the motion of that parcel and the chemical evolution is only apparent from semi-analytical 1D disk models. Because of the algebraic nature of equations for the disk structure, the kinetics can be attached quite easily to the parcel of gas within the accretion flow as has been done by Duschl et al. (1996), Bauer et al. (1997), Finocchi et al. (1997), Finocchi & Gail (1997), and Aikawa et al. (1999). However, with respect to the numerical scheme, one should exercise caution when integrating the system of stiff differential equations which is mainly done by implicite multi-step techniques such as the BDF method (cf. Strehmel & Weiner 1995). In addition to the stability of the numerical scheme, implicite multi-step techniques take advantage of increasing step sizes in time, which are in general many orders of magnitude longer than the initial step size. Considering only one single integration step, the shift between the previous point t_n and the new point t_n+1 can become quite large. This means that important physical processes, taking place meanwhile, may not be resolved because of the numerically determined step size control and order selection.

Because of the non-algebraic nature of equations, determining the 1+1-D-disk structure, the chemical and the physical disk evolution along the streamline cannot be solved simultaneously as it was done in the case of semianalytical 1D-disk models. Therefore, some kind of reasonable decoupling of kinetics and advective transport is needed.

Returning to the extended discussion concerning the implementation of a Lagrangian fluid element in a 1D-chemical model, we are able to specify a condition for decoupling quite easily. Adopting a cylindrical coordinate system, the spiral motion of the Lagrangian fluid element towards the central object can be analysed in terms of the vectorial components of the accretion flow. Concerning the azimuthal component, the motion is assumed to be on Keplerian orbits. It can be shown that the azimuthal component of the velocity is highly supersonic while the value of the radial component is much less than the sound speed. Therefore, it seems reasonable to decompose the trajectory into a Keplerian orbit which is possessed of a small radial drift. In the case of steady 1+1-D-$\alpha$-disk models the radial component of the velocity is always directed inwards.

Following the trajectory within two coaxial cylinders with radii r₁, r₂ we impose the condition that the variation of the radial component of the accretion flow $v_{\rm r}$ vanishes vertically. This assumption is in agreement with the thin-disk approximation determining $v_{\rm r}$ by Eq. (2). Because of the assumed Keplerian motion, we have to take into account only the radial component of the accretion flow. As long as the physical conditions within the two radii are quite similiar there is no need to sustain the coupling of kinetics and advective transport. In particular, those physical quantities important in influencing the chemical evolution are the temperature T(r,z), the mass density $\varrho(r,z)$, and the optical depth $\tau(r,z)$. The decoupling mechanism was applied on a distance $\Delta r\!=\!r_2-r_1$ where the local change within the T, $\varrho$, and $\tau$ profiles is less than 10%. For example, the decomposed trajectory of the simulation mod1ELIS is shown in Fig. 2.

	Figure 2: Typical decomposition of a trajectory within a steady 1+1-D-$\alpha$-disk model from 10 to 1 AU; mod1ELIS. More than 366 gridpoints are needed to construct the grid. The change in the physical quantities moving from one gridpoint to the neighboring gridpoint is 9% on average.
Open with DEXTER

After calculating the time t needed to follow the transport from r₁ to r₂, we first solved the kinetic equations within the vertical domain at r₁ for $\Delta t\!=\!t$. Second, the fluid element was transported by advection to r₂ taking the equation of continuity into account.

We calculated the chemical evolution along three different trajectories from 10 AU to 1 AU with the parameters given in Table 1. In the last plot of Fig. 3, the time needed for advective transport is plotted vs. the radial position. As we shall see, the results presented in the next section can be classified according to local disk structure and trajectory. With respect to the chemical evolution comparisons can be made of i) similiar local disk structure and different trajectories and ii) different local disk structure and similiar trajectories. The latter case refers to simulations mod2ELIS and mod3ELIS while mod1ELIS and mod3ELIS represent the former.

Table 1: Parameter of the presented simulations determining the viscous evolution of the disk.

4.2 Diffusion

The previous discussion might lead one to separate the kinetics and vertical transport processes. However, there is a subtle difference concerning the transport properties. In a multi-component system, the various species will move due to diffusion at different diffusion velocities, cf. Eqs. (13) and (14). In particular, the bulk velocity of the mixture is determined by the local composition of the mixture itself. Splitting techniques such as operator splitting were applied, e.g., by Rawlings & Hartquist (1997), in order to solve the reaction-diffusion model. However, the problems related to splitting strategies are apparent from the lack of an uniform time scale, independent of the kinetics.

Instead of decoupling of kinetics and diffusion one can integrate Eq. (16) simultaneously. The discretised equations are given by semi-discretisation of the vertical components of Eq. (16) adapted for a Lagrangian description. Finally, the system of n partial differential equations was transformed into a set of $n \times nz$ ordinary differential equations where n and nz denote the number of species and the number of vertical gridpoints,

	Figure 3: Temperature distribution of a steady 1+1-D-$\alpha$-disk model refering to 3 different simulations mod1ELIS, mod2ELIS, and mod3ELIS (clockwise from top left). The temperature contours for 100 K and 500 K are shown. The figure on the right bottom refers to the trajectories of each simulation within the computational domain.
Open with DEXTER

respectively. Although the dimension of the system may easily become quite large, the system itself is characterized by a sparse matrix. Depending on the used discretisation scheme for approximating the spatial derivatives, the system is determined by a banding structure at least.

Because of the global transport processes within a reaction-diffusion model, boundary conditions have to be fixed. Apparently the condition refering to the disk's midplane ($z\!=\!0$) is determined by the geometry of the disk itself. Furthermore, we assumed that the gradients in the species' abundances vanish beyond the photospheric height $z_{}^{\ast}$.

4.3 Chemical network

The disk chemical model is based on that of Markwick et al. (2002). The species set consists of a total of $n\!=\!242$ species including 69 mantle components made up of the elements H, ${\rm He}$, C, N, O, S, ${\rm Si}$, ${\rm Fe}$ and ${\rm Mg}$. About 2500 reactions connect these species.

The model contains various types of reactions in the gas phase: ion-neutral and neutral-neutral reactions as well as photoprocesses. The interaction between the species in the gas phase and the chemically reactive mantle components onto grain surfaces is described by adsorption and desorption mechanisms. Adsorption is caused by freezing of species onto grains. The reverse processes, thermal desorption and cosmic ray heating, desorb species from grain surfaces.

Expressions for the rate coefficients $k_{}^{\pm}$, refering to reactions in the gas phase as well as to adsorption processes onto grains, can be obtained by macroscopic considerations. If $\sigma_{{\rm c}}$ and $\vec{v}$ denote the cross-section and the velocity of the molecules, respectively, the mean collision rate coefficient is given by $\!<\!\!\sigma_{{\rm c}} \vec{v}\!\!>\!$. Once averaged over a Maxwellian velocity distribution, one gets the Arrhenius type of rate coefficients which is typically used in astrochemical databases. For the adsorption process, one has additionally to take the sticking probability into account which depends on the charge of the adsorbing particle, cf. Umebayashi & Nakano (1980).

The mobility of the adsorbed species on the surface depends on the surface structure. Therefore, in addition to bimolecular processes one has to consider other microscopic effects such as the characteristic (lattice)-vibration frequency in order to describe the rate coefficients related to desorption processes on grains. With respect to thermal desorption as well as to cosmic ray heating the kinetic equations are a variation of the Polanyi-Wigner equation, cf. Hasegawa et al. (1992) and Hasegawa & Herbst (1993).

Sources of ionisation in the disk include energetic particles formed by the radioactive decay of $_{}^{26}{\rm Al}$, interstellar UV photons, and ionisation due to cosmic rays. The rate of the latter process depends on the local value of the surface density $\Sigma(r,z)$, cf. Umebayashi & Nakano (1981).

5 Results

The chemical evolution was evaluated along the trajectory from 10 ${\rm AU}$ to 1 ${\rm AU}$. The decomposition of the trajectory needed to decouple kinetics and advective transport requires more than 366 gridpoints radially (mod1ELIS, cf. Fig. 2). The number of gridpoints in the other simulations mod2ELIS and mod3ELIS have the same order of magnitude. The computational domain was divided vertically into an equidistant grid with $nz\!=\!40$. As

	Figure 4: Chemical evolution of ${\rm CS}$ and ${\rm SO_2}$ along the trajectory from 10 AU to 1 AU; mod1ELIS without vertical mixing. The contour lines refer to the abundance values relative to the total amount of hydrogen: 10-10,-13,-16 and additionally 10-7 for ${\rm SO_2}$.
Open with DEXTER

soon as vertical mixing was included, the Crank-Nicholsen scheme was applied to discretise the spatial derivatives. Finally, the vode-integrator developed by Brown et al. (1988) was used to solve the system of nonlinear ordinary differential equations. The initial abundances of the species at 10 ${\rm AU}$ were taken from Willacy et al. (1998).

The presentation of our results will be made on three different levels in order to highlight different aspects of the chemical evolution in disks. On the first level, we will analyse the data obtained by a simulation without vertical mixing. To some extent we will use this to identify the main formation and destruction processes of particular species. Once the chemical evolution is determined by a reaction-diffusion model, it is easier to draw conclusions concerning the influence of mixing processes on the chemical evolution (second level). Finally, we will clarify the influence of the disk dynamics on the chemical evolution. On this third level, we analyse the chemical evolution with respect to different disk parameters determining the dynamics of mod1ELIS-mod3ELIS, while the disk parameters have been kept constant for the previous discussion.

At each level, we will consider the chemical evolution of the ${\rm CS}$ molecule. Because of its large dipole moment, ${\rm CS}$ is expected to trace high-density regions. Based on measurements of higher rotational ${\rm CS}$ transitions one might get better constraints on the chemical and physical structure of disks. In addition, since ${\rm CS}$ is detected in the atmospheres of comets - cf., e.g., Snyder et al. (2001) - one may be able to compare the early evolutionary stages of protostellar disks and the early stages of the Solar Nebula. The ongoing use of ${\rm CS}$ as a diagnostic of different star-forming regions should be seen against this background - cf., e.g., Blake et al. (1992), Dutrey et al. (1997), Launhardt et al. (1998), Qi (2000), and Doty et al. (2002).

The analysis of ${\rm CS}$ cannot be decoupled from the chemistry of other sulphur-bearing molecules. Therefore, the evolution of other species closely connected to the evolution of ${\rm CS}$ such as ${\rm SO}$, ${\rm SO_2}$, ${\rm HCS_{}^+}$, ${\rm HCO_{}^+}$, and ${\rm NH_3}$ has to be discussed as well.

5.1 Chemical evolution based on a pure kinetic model

We start the discussion with the question of where most of the sulphur is locked. Sulphur-bearing molecules may stick on the dust surfaces due to adsorption processes. In particular, sulphur is mainly channSte eled into mantle components in low temperatures regions. Note that in general more than 50% of the reservoir of sulphur is assumed to be bound within the grain core itself, especially in Troilite (${\rm FeS}$), cf. Pollack et al. (1994) and Finocchi et al. (1997). For increasing temperatures, one obtains a rapid rise of sulphur-bearing molecules in the gas phase. Now most of the sulphur is locked in ${\rm CS}$ and ${\rm SO_2}$ molecules, cf. Fig. 4.

The conclusions drawn by Doty et al. (2002), Willacy et al. (1998), and Langer et al. (2000) are confirmed by our calculations. Doty et al. emphasise that sulphur is mainly channeled into ${\rm SO_2}$ in regions with $T\!\approx\! 100$- $400 \ {\rm K}$ while ${\rm CS}$ serves as a good measure of sulphur in the cool exterior regions. The simulations of Willacy et al., which take only the evolution along the disk midplane into account, have shown that ${\rm SO_2}$ is the dominant sulphur-bearing species. In addition, because of the Lagrangian description, we may characterise the high abundance of ${\rm CS}$ by a disk chronology analogy (Langer et al. 2000). They argued that an increased abundance of ${\rm CS}$ is associated with early states of chemical evolution.

We will now follow the trajectory from 10 to 1 AU, which crosses regions with different physical conditions. At first the relative abundance of ${\rm CS}$ decreases while following the accretion flow toward inner regions. In particular, in the colder temperature regions, close to the photospheric disk height, the amount of available sulphur in the gas phase drops due to adsorption onto the grain surfaces. For the vertically adjacent regions, the temperature will increase. Thus, the efficiency of reactions in the gas phase rises. Compared to the initial abundance at $10 \ {\rm AU}$ carbon monosulfide becomes more abundant in regions represented by the contour lines of 10^-16 and 10^-13, cf. Fig. 4. The major routes to ${\rm CS}$ in the disk are related to the evaporation of ${\rm H_2CS}$ and ${\rm H_2S}$ from grain mantles. Reaction of these molecules with ${\rm C_{}^+}$ form ${\rm CS}$ and ${\rm HCS_{}^+}$, respectively. In addition, proton transfer reactions with ${\rm H_2CS}$ followed by the dissociative recombination of ${\rm H_3CS_{}^+}$ also produces ${\rm CS}$. In the inner region of the disk, the most important reactions involve ${\rm HCS_{}^+}$

$$HCS_{}^+ + e_{}^- \longrightarrow CS + H.$$ (R 1)

Because of the inverse temperature dependence of its rate coefficient ( $k\!\propto\!T_{}^{-3/4}$), reaction (R 1) gives rise to higher abundances of ${\rm CS}$ especially in colder regions. Of course, an enhanced abundance of ${\rm HCS_{}^+}$ needed to produce ${\rm CS}$ can be found in these regions, cf. Fig. 5. The destruction of ${\rm CS}$ via protonation of ${\rm HCO_{}^+}$

$$HCO_{}^+ + CS \longrightarrow HCS_{}^+ + CO$$ (R 2)

cannot compensate the production of ${\rm CS}$ via reaction (R 1).

Once one passes into regions with temperatures above 100 K, desorption processes become more efficient and an enriched gas phase chemistry can be expected. An increased abundance of ${\rm CS}$ within the contour line of 10^-10 (cf. Fig. 4) is attributed to the evaporation of mantle molecules back to the gas phase. Of course, the adsorbed counterpart of ${\rm CS}$, ${\rm CS_{{\rm ads.}}}$, becomes less abundant, cf. Fig. 5.

In addition to desorption, the change in the ${\rm CS}$ abundance can be explained by differences within the gas phase chemistry. The destruction of ${\rm CS}$ is influenced by oxygen chemistry. Because of the enhanced abundance of atomic oxygen (cf. Fig. 5, region between 3 and 5 AU) ${\rm CS}$ is efficiently removed by

$$CS + O \longrightarrow CO + S$$ (R 3)

returning sulphur to atomic form. Similar dependencies were reported by Charnley (1997) and Doty et al. (2002), although different environments were considered. However, most of the available sulphur in the gas phase is now locked in ${\rm SO_2}$. The increased abundance of ${\rm CS}$ for $r\!\le\!2$ AU is related to the enhanced abundances of ammonia, ${\rm NH_3}$, which has a very large proton affinity. The proton transfer reaction between ammonia and the thioformyl cation, ${\rm HCS_{}^+}$,

$$NH_3 + HCS_{}^+ \longrightarrow CS + NH_4^+$$ (R 4)

reforms ${\rm CS}$. Note that the corresponding rate coefficient is independent of the temperature. It is remarkable that ammonia can contribute

Figure 5: Chemical evolution of ${\rm HCS_{}^+}$, ${\rm O}$, ${\rm S}$, ${\rm NH_3}$, ${\rm CS_{{\rm ads.}}}$, and ${\rm HCO_{}^+}$ (clockwise from top left) along the trajectory from 10 to 1 ${\rm AU}$; mod1ELIS without vertical mixing. The contour lines refer to abundance values relative to the total amount of hydrogen of 10-20,-21,-22 ( ${\rm HCS_{}^+}$), 10-13.5 (${\rm O}$), 10-15,-16 ( ${\rm HCO_{}^+}$), 10-15,-17 (${\rm S}$), 10-15,-16,-18 ( ${\rm CS_{ads}}$), and 10-5,-6 ( ${\rm NH_3}$).

Open with DEXTER

significantly to the increased abundance of ${\rm CS}$ although the coupling between nitrogen and sulphur chemistries, through ${\rm NS}$, is weak. The implications of this may be of particular interest to the unexpected distribution of ${\rm CS}$ and ${\rm NH_3}$ seen in the observations of a few star-forming regions, cf. Taylor et al. (1996).

5.2 Chemical evolution based on a reaction-diffusion model

The previous model demonstrated that vertical abundance gradients of the gas phase species develop. Therefore, the gradients of the chemical potentials cause a diffusion flow.

	Figure 6: Chemical evolution of ${\rm CS}$, ${\rm SO}$, ${\rm SO_2}$, and ${\rm S}$ (clockwise from top left) along the trajectory from 10 to 1 ${\rm AU}$; mod1ELIS incl. vertical mixing. The contour lines refer to abundance values relative to the total amount of hydrogen of 10-13,-16 and 10-10 for ${\rm SO_2}$ additionally.
Open with DEXTER

The simulations in this subsection differ from the previous simulations only in respect to the inclusion of the reaction-diffusion process. Neither the chemical network including the initial abundances nor the disk parameters were changed.

At the beginning, the chemical evolution along the trajectory is quite similiar to that with vertical mixing neglected. In the particular case of ${\rm CS}$, differences in the vertical distribution become significant only for $r \! \le \! 8.5 \ {\rm AU}$, cf. Fig. 6. For $r \! > \! 8.5 \ {\rm AU}$ diffusion processes cannot really influence the chemical evolution significantly because of the low abundances of molecules in the gas phase. Nevertheless, the already established gradients are removed by diffusion. Apart from the gradients, the diffusion velocity also depends on the local value of the diffusion coefficient. Concerning ${\rm CS}$, one expects to observe the most effective mixing near the disk midplane, thus transporting ${\rm CS}$-molecules to regions above the disk midplane. Here, because thermal adsorption processes are favoured at

	Figure 7: Chemical evolution of ${\rm CS_{ads}}$, ${\rm O}$, ${\rm HCS_{}^+}$, and ${\rm NS}$ (clockwise from top left) along the trajectory from 10 to 1 ${\rm AU}$; mod1ELIS incl. vertical mixing. The contour lines refer to abundance values relative to the total amount of hydrogen of 10-15,-16 ( ${\rm CS_{ads}}$), 10-13.5 (${\rm O}$), 10-22.5 (${\rm NS}$), and 10-20,-21,-22 ( ${\rm HCS_{}^+}$).
Open with DEXTER

low temperatures, the molecules will freeze out onto grain surfaces. Therefore, the contour lines of the adsorbed counterpart ${\rm CS_{ads.}}$ are shifted to higher values of z, cf. Fig. 6. There is no enrichment of ${\rm CS}$ between 8 and $9 \ {\rm AU}$ as was found in the pure kinetic model, cf. Fig. 5. In comparision to the results of the previous level, less ${\rm CS}$ molecules were adsorbed onto grains in these regions, and thus less mantle material is available for desorption.

Between 4 and $6 \ {\rm AU}$ the destruction of ${\rm CS}$ in the gas phase via reaction (R 3) becomes more effective because of the increased abundance of atomic oxygen, cf. Fig. 7. The oxygen itself is formed mainly by the reaction

$$N + NO \longrightarrow N_2 + O.$$ (R 5)

Unlike the previous model, the destruction of ${\rm CS}$ takes place in the inner disk regions although there is ammonia available at $r \! \le \! 1.5 \ {\rm AU}$ to recycle ${\rm CS}$ via reaction (R 4). However, because of the low ${\rm HCS_{}^+}$ abundance, cf. Fig. 7, the formation rates are not sufficient to overcome the loss rates with atomic oxygen. With respect to a chronological classification one may argue, due to the earlier appearance of the peak in ${\rm CS}$, that ${\rm CS}$ and ${\rm NH_3}$ are associated with earlier and later stages of the chemical evolution, respectively. Langer et al. (2000) have drawn the same conclusions.

Finally, relating to the previous model, we discuss which sulphur-bearing species in the gas phase contributes most to the entire reservoir of sulphur. Compared to the results of the previous model, cf. Fig. 4, ${\rm SO_2}$ seems to be already a good measure of sulphur abundance in the gas phase at larger radii. Additionally, the abundance of ${\rm SO_2}$ molecules drops in the very inner regions in sharp contrast to the previous model because of the low abundance of the ${\rm OH}$ radical which is more easily forced into ${\rm H_2O}$ than into ${\rm SO_2}$ through the reaction

$$SO + OH \longrightarrow SO_2 + H.$$ (R 6)

On the other hand, the abundance of atomic sulphur does not decrease due to the contributions of the reactions (R 3) and

$$N + SO \longrightarrow NO + S$$ (R 7)

$$N + NS \longrightarrow N_2 + S.$$ (R 8)

5.3 Chemical evolution along different trajectories

The previous results clearly demonstrate mixing processes do influence the chemical evolution. Of course, the local change in the composition of species due to diffusion may have an impact on the chemical evolution globally since the disk material is transported by advection. Given that we have a reaction-diffusion model to determine the kinetics, we now focus on how the advective mass transport acts on the chemical evolution of the disk.

Because of the steady-state assumption, the dynamics of a protoplanetary accretion disk are represented only by $\alpha$ and $\dot{M}$. Therefore, discussions of the disk dynamics become only reasonable if we compare different steady-state solutions in terms of the parameters $\dot{M}$ and $\alpha$. The trajectories as well as the local disk structures referring to three different steady-state solutions are shown in Fig. 3.

We now compare the models mod1ELIS and mod3ELIS (see Fig. 8). They can be classified by different trajectories, but quite similar disk structures. We find no significant differences in the chemical evolution of ${\rm CS}$. Although the vertical distribution of ${\rm CS}$ is slightly larger in mod3ELIS than in mod1ELIS, this reflects the fact that the vertical extent of the former disk is larger than that of the latter. In particular, the main formation and destruction processes along the trajectories are identical. However, one obtains a gradual shift in the abundances of, e.g., ${\rm SO}$, ${\rm SO_2}$, and ${\rm S}$. Although the evolution of these molecules follows the same chronological sequence in both simulations, there is still a change in the abundances in relative high temperature regions as well as in regions with lower temperatures. The change of the position within the contour lines does not follow the physical conditions such as, e.g., the temperature profile. Therefore, different aspects have to be taken into account such as, e.g., minor differences in the temperature structure as well as different transport velocities in general.

Let us consider, for instance, the contour line of ${\rm SO}$ referring to the value of 10^-13 relative to the total amount of hydrogen. Compared to the result of mod1ELIS, the contour line of mod3ELIS is shifted to larger radii (Figs. 6 and 8). The local variation in the abundances can be best explained

	Figure 8: Chemical evolution of ${\rm CS}$, ${\rm SO}$, ${\rm S}$, and ${\rm SO_2}$ (clockwise from top left) along the trajectory from 10 AU to 1 AU; mod3ELIS including vertical mixing. The contour lines refer to the abundance values relative to the total amount of hydrogen: 10-13,-16 and additionally 10-10 for  ${\rm SO_2}$.
Open with DEXTER

by time dependent chemical effects rather than by minor differences in the efficiency of adsorption and desorption processes due to slightly different temperature profiles. That is, the time needed for transport between $10 \ {\rm AU}$ and $8.5 \ {\rm AU}$ along the trajectory of mod1ELIS is about $20~000 \ {\rm yrs}$ while the radial transport of mod3ELIS operates on a timescale roughly twice as long, cf. Fig. 3. Therefore, in the case of mod1ELIS less ${\rm SO}$ mantle material can be ejected into the gas phase.

In the outer regions of the computational domain, a small variation within the temperature profiles can only account for differences in the abundances if the thermal desorption is affected sensitively by local changes in the temperature. This, of course, depends on the species-dependent value of binding energies, and therefore species with high binding energies, such as ${\rm SO_2}$, are expected to show temperature-dependent variations in their abundances in those regions. Additionally, in terms of the ${\rm SO_2}$ sublimation temperature of $90 \ {\rm K}$ (Wutz et al. 1992) the shift within the contour lines, e.g., of 10^-10, cf. Figs. 6 and 8, respectively, follows exactly the variation of the temperature profile.

It was already pointed out that the chemical evolution based on the simulations mod1ELIS and mod3ELIS can be described by the same sequence of main formation and destruction processes. However, important differences within the relative abundances can still be found,

	Figure 9: Chemical evolution of ${\rm OH}$, ${\rm O}$, and ${\rm O_2}$ along the trajectory within 4 AU; The figures on the left and right panel refer to results of the simulations mod1ELIS and mod3ELIS, respectively. Vertical mixing were taken into account in both simulations.
Open with DEXTER

especially in regions with $r \! \le \! 5 \ {\rm AU}$. The different values in the abundances mainly originate from different timescales which the radial transport operates on. As long as one compares the chemical evolution within the same temper ature region, the evolution is quite similiar. For example, the destruction of ${\rm SO}$ which occurs in both simulations for $T\!\ge\!300\ {\rm K}$ is illustrative of this particular behaviour.

The chemical evolution of ${\rm SO}$ is intimately connected to the oxygen chemistry and, in particular, with the distribution of reactive oxygen. The "destruction'' of ${\rm SO}$ for $r \! \le \! 3 \ {\rm AU}$ is mainly caused by missing production routes such as

$$S + OH \longrightarrow SO + H$$ (R 9)

and

$$S + O_2 \longrightarrow SO + O.$$ (R 10)

While the abundance of atomic sulphur remains approximately constant, most of the hydroxyl will be channeled into water via

$$OH + H_2 \longrightarrow H_2O + H$$ (R 11)

at higher temperatures. For the temperature region $T\!\ge\!300\ {\rm K}$, the formation of ${\rm OH}$ is mainly due to the reaction

$$O + H_2 \longrightarrow OH + H.$$ (R 12)

Consequently, significant production of molecular oxygen in the gas phase via the reaction

$$O + OH \longrightarrow O_2 + H.$$ (R 13)

is inhibited since the amount of ${\rm O}$ and ${\rm OH}$ is

	Figure 10: Chemical evolution of ${\rm CS}$, ${\rm SO}$, ${\rm S}$, and ${\rm SO_2}$ (clockwise from top left) along the trajectory from 10 AU to 1 AU; mod2ELIS including vertical mixing. The contour lines refer to the abundance values relative to the total amount of hydrogen: 10-13,-16 and additionally 10-10 for ${\rm SO_2}$.
Open with DEXTER

limited at increasingly higher temperatures.

Although the kinetics interacts with advective and diffusive mass transport processes the chemical evolution may be governed by local properties like the temperature structure. The explanation for the destruction of ${\rm SO}$ given above holds for both simulations mod1ELIS and mod3ELIS. The quite similiar distributions of ${\rm O}$, ${\rm O_2}$, and ${\rm OH}$ along the trajectories, cf. Fig. 9, may be explained in terms of the same chemical signature determined mainly by the temperature structure. Therefore, in addition to the importance of global transport processes in disks on the chemical evolution, we should never underestimate the importance of the parameters of the kinetics.

Continuing the discussion, we now refer to the simulations mod2ELIS and mod3ELIS which have quite different local disk structure but with approximately similar transport by advection (cf. Figs. 8 and 10). Because of the minor differences between the time scales of advective transport below $t \! \le \! 70~000 \ {\rm yrs}$, cf. Fig. 3, the differences within the local disk structures of mod2ELIS and mod3ELIS beyond $r \! \ge \! 7 \ {\rm AU}$ account for the different chemical evolution. This region is of special interest as reactions in the gas phase begin to dominate the chemical evolution in mod3ELIS while the gas-grain interactions still remain the driving source in mod2ELIS. However, peaks in the abundances of species in the gas phase, signifying the contributions due to thermal desorption processes, are linked to temperature regions where these processes become effective.

It is remarkable that with respect to the late stages, that is in the inner disk, the chemical evolution based on mod2ELIS and mod3ELIS, respectively, is quite similar and no additional reactions have to be taken into account to explain the chemical evolution.

6 Conclusion

We have investigated the chemical evolution of a protoplanetary disk. For this study, we used the 1+1-D-approximation of cylindrical steady disk models. Using a Lagrangian description, the change in the molecular composition was evaluated along different trajectories. Besides advective transport we also considered transport processes by diffusion.

As these models have been prepared for the study of a reacting flow system, we analysed the underlying physical transport properties in detail in order to estimate where coupling of mass transport processes and kinetics is needed. Modeling the vertical mixing requires a reaction-diffusion model while in the case of the advective transport decoupling seems to be an appropriate approach to treat both processes.

The global chemical evolution was found to change in all models due to vertical mixing processes, compared to results without vertical mixing. Furthermore, the local changes in abundances due to diffusion will be "transported'' by advection. Diffusion processes appear to have only a limited effect on the chemical evolution in regions where gas-grain proccesses dominate the kinetics.

In particular, models with different disks parameters $\dot{M}$ and $\alpha$ may show similarities in the chemical evolution of the disk gas. Thus, in addition to mass transport processes, information about the local conditions, which determine the kinetics, are still needed in order to analyze the chemical evolution.

Appendix A: Boundary conditions determined by radiation transport

For a one-dimensional planar atmosphere, the steady radiation transfer equation is given by

 

$$\cos\Theta \frac{\partial I_{\nu}}{\partial \tau_{\nu}} I_{\nu}(\tau_{\nu},\Theta)-S_{\nu}(\tau_{\nu})$$ = (A.1)

$$\tau_{\nu} \int\limits_{-\infty}^{z} \kappa_{\nu}(z)\ {\rm d}z,$$ =

where $\Theta$, $I_{\nu}$, $S_{\nu}$, $\tau_{\nu}$, and $\kappa_{\nu}$ refer to the photon propagation angle and the frequency-dependent quantities: the specific intensity of radiation, the source function, the optical depth, and the absorption coefficient. z denotes the geometrical depth. The emergent specific intensity in general and its value outside the atmosphere, i.e., $\tau_{\nu}\!=\!0$, can be obtained by

 

$$I_{\nu}(\tau_{\nu},\Theta) \sec \Theta \int\limits_{0}^{\infty} S_{\nu} (\tau_{\nu}^{\ast}) ... ...}_{}^{ -(\tau_{\nu}^{\ast}-\tau_{\nu}) \sec \Theta } \ {\rm d}\tau_{\nu}^{\ast}$$ = (A.2)

$$I_{\nu}(0,\Theta) \sec \Theta \int\limits_{0}^{\infty} S_{\nu} (\tau_{\nu}^{\ast}) {\rm e}_{}^{ - \tau_{\nu}^{\ast} \sec \Theta } \ {\rm d}\tau_{\nu}^{\ast} .$$ =

In particular, the monochromatic radiation flux $\pi F_{\nu}$ can be specified by

 

$$F_{\nu}(\tau_{\nu}) 2\int\limits_{0}^{\pi} I_{\nu}(\tau_{\nu},\Theta) \cos \Theta \sin \Theta \ {\rm d}\Theta$$ = (A.3)

$$F_{\nu}(0) 2\int\limits_{0}^{\infty} S_{\nu} (\tau_{\nu}^{\ast}) K_2(\tau_{\nu}^{\ast}) \ {\rm d}\tau_{\nu}^{\ast},$$ =

where the last equation corresponds to the flux outside the atmosphere. K₂ denotes the exponential integral of the second kind which belongs to the family

$$K_n(x)=\int\limits_1^{\infty} y_{}^{-n} {\rm e}_{}^{-xy} {\rm d}y.$$

Replacing the source function $S_{\nu}$ by its power expansion i.e.,

$$S_{\nu}(x) = S_{\nu}(x_{0}) \left. \frac{\partial S_{\nu}}{\parti... ... \right\vert _{x_{0}} \! \! \! (x-x_{0})_{}^2 + \mathcal{O}\left(x_{}^3\right),$$

substituting into Eq. (A.3), and taking the following expression into account

$$\int\limits_{0}^{\infty}x_{}^s K_n(x) \ {\rm d}x = \frac{s!}{s+n},$$

one obtains the final expression for the monochromatic radiation flux $\pi F_{\nu}(0)$ (Unsöld 1968):

 

$$F_{\nu}(0) = S_{\nu}(\tau_{\nu}^{\ast})+\left. \left( \frac{2}{3}... ...{\nu}^2} \right\vert _{\tau_{\nu}^\ast} + \mathcal{O}\left(\tau_{\nu}^3\right).$$ (A.4)

Neglecting the nonlinear terms in (A.4), we find $F_{\nu}(0)~=~S_{\nu}\ (\tau_{\nu}^{\ast}~=~2/3)$, which is well-known as the Eddington-Barbier relation. The emergent monochromatic radiation flux is equal to the source function at optical depth $\tau_{\nu}^{\ast}\!=\!2/3$.

With respect to the vertical disk model, the transport of radiation energy has to be modified within the frequency-integrated variables I, S, and $\tau$. With exception of grey material, i.e., $\kappa\equiv\kappa_{\nu}$, the transfer equation obtained serves only as an approximation. Providing a frequency-averaged absorption coefficient in terms of a single mean opacity $\overline{\kappa}\!=\!\kappa$, one obtains the zeroth and first angular moment of the transfer Eq. (A.1):

  

$$\frac{1}{\kappa_{(\nu)}} \frac{\partial \ }{\partial z} H_{(\nu)} J_{(\nu)} - S_{(\nu)}$$ = (A.5)

$$\frac{1}{\kappa_{(\nu)}} \frac{\partial \ }{\partial z} K_{(\nu)} H_{(\nu)},$$ = (A.6)

determinded by the frequency-integrated variables S, J, H, and K. The angular moments of the (monochromatic) specific intensity are given by

$$\left( J_{\nu}, H_{\nu}, K_{\nu} \right) = \int I_{\nu} \lef... ...{}^2 \Theta \right) \frac{{\rm d}\omega}{4\pi}\cdot \nonumber$$

In addition, the indices within the parentheses refer to the exact (frequency-dependent) description. Based on this set of equations one can derive a single mean opacity $\overline{\kappa}$ which guarantees the correct total flux of radiation energy, i.e. $\int H_{\nu} {\rm d}\nu = H$,

$$\int\limits_{0}^{\infty} \frac{1}{\kappa_{\nu}} \frac{\parti... ..._{0}^{\infty} \frac{\partial K_{\nu}}{\partial z} {\rm d}\nu .$$ (A.7)

The assumption that the specific intensity $I_{\nu}$ is close to its blackbody value holds as long as one considers regions of a non-negligible amount of optical depth. Because of the isotropic radiation field, the $\cos_{}^2 \Theta$ term can be replaced by its mean value 1/3 which finally results in Eddington's approximation $K_{(\nu)}=1/3 \ J_{(\nu)}$ and the Rosseland mean opacity $\overline{\kappa}_{{\rm R}}$, respectively. This kind of mean opacity was used to approximate the radiative energy transport in Eq. (7). Corresponding to the relationship (A.7) the boundary condition for $\tau_{\nu}^{\ast}\!=\!2/3$ can be replaced by $\tau_{}^{\ast} \!=\! 2/3$.

Finally, we have to specify the relation between the physical temperature T and the effective temperature $T_{{\rm eff}}$ used to measure the emitted flux from a black body, see Eq. (5). By requiring radiative equilibrium, Eqs. (A.5-A.6) yield (Unsöld 1968)

$$J=3H\left( \frac{2}{3}+\tau \right),$$ (A.8)

where Eddington's approximation was used. Applying the Stefan-Boltzmann law one obtains the simple relation

$$T_{}^4=\frac{3}{4}\left( \frac{2}{3}+\tau \right) T_{{\rm eff}}^4 \ .$$ (A.9)

With respect to the chosen boundary value of $\tau_{}^{\ast}$ one can easily specify $T(z_{}^{\ast})\!=\!T_{{\rm eff}}$.

Acknowledgements

M.I. thanks Edwin Bergin for useful discussions. We thank P. D'Alessio for providing data of her models. M.I. was supported by DFG grant He 1935/17-1. Research in Astrophysics at UMIST is supported by PPARC. The computations were performed, using the equipment of the Department of Mathematics, Friedrich-Schiller-Universität Jena.

References

 
• Aikawa, Y., Umebayashi, T., Nakano, T., & Miyama, S. 1999, ApJ, 519, 705 In the text NASA ADS  
• Aikawa, Y., & Herbst, E. 1999, A&A, 351, 233 In the text NASA ADS  
• Aikawa, Y., & Herbst, E. 2001, A&A, 317, 1107 In the text  
• Aikawa, Y., van Zadelhoff, G., van Dishoeck, E., & Herbst, E. 2002, A&A, 386, 622 In the text NASA ADS  
• Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 In the text NASA ADS  
• Bauer, I., Finocchi, F., Duschl, W., Gail, H.-P., & Schlöder, J. 1997, A&A, 317, 273 In the text NASA ADS  
• Bell, K. R., Cassen, P., Klahr, H. H., & Henning, Th. 1997, ApJ, 486, 372 In the text NASA ADS  
• Blake, G., van Dishoeck, E., & Sargent, A. 1992, ApJ, 391, L99 In the text NASA ADS  
• Brown, P., Byrne, G., & Hindmarsh, A. 1988, LLNL Report UCRL-98412 In the text  
• Charnley, S. 1997, ApJ, 481, 396 In the text NASA ADS  
• Chiang, E., & Goldreich, P. 1997, ApJ, 490, 368 In the text NASA ADS  
• D'Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411 NASA ADS  
• Doty, S., van Dishoeck, E., van der Tak, F., & Boonman, A. 2002, A&A, 389, 446 In the text NASA ADS  
• Duschl, W., Gail, H.-P., & Tscharnuter, W. 1996, A&A, 312, 624 In the text NASA ADS  
• Dutrey, A., Guilloteau, S., & Guélin, M. 1997, A&A, 317, L55 In the text NASA ADS  
• Finocchi, F., Gail, H.-P., & Duschl, W. 1997, A&A, 325, 1264 In the text NASA ADS  
• Finocchi, F., & Gail, H.-P. 1997, A&A, 327, 825 In the text NASA ADS  
• Gail, H.-P. 1998, A&A, 332, 1099 In the text NASA ADS  
• Gail, H.-P. 2001, A&A, 378, 192 In the text NASA ADS  
• Gail, H.-P. 2002, A&A, 390, 253 In the text NASA ADS  
• Hasegawa, T., Herbst, E., & Leung, C. 1992, ApJ, 82, 167 In the text  
• Hasegawa, T., & Herbst, E. 1993, MNRAS, 263, 589 In the text NASA ADS  
• ayashi, C., Kakazawa, N., & Nakagawa, Y. 1985 in Protostars and Planets II , ed. D. Black & M. Matthews (Arizona Press), 1100 In the text  
• Kley, W., & Lin, D. N. C. 1992, ApJ, 397, 600 In the text NASA ADS  
• Kley, W., Papaloizou, J. C. B., & Lin, D. N. C. 1993, ApJ, 416, 679 In the text NASA ADS  
• Langer, W., van Dishoeck, E., Bergin, E., et al. 2000, in Protostars and Planets IV, ed. V. Mannings & A. Boss, 29-58 (Tucson: University of Arizona) In the text  
• Launhardt, R., Evans, N., Wang, Y., et al. 1998, ApJS, 119, 59 In the text NASA ADS  
• Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 In the text NASA ADS  
• Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 In the text NASA ADS  
• Markwick, A., Ilgner, M., Millar, T, & Henning, Th. 2002, A&A, 385, 632 In the text NASA ADS  
• Pollack, J. B., Hollenbach, D., Beckwith, S., et al. 1994, ApJ, 421, 615 NASA ADS  
• Pringle, J. E. 1981, A&ARv, 19, 137 In the text  
• Qi, C. 2000, Ph.D. Thesis, California Institute of Technology, Pasadena In the text  
• Rawlings, J., & Hartquist, T. 1997, ApJ, 487, 672 In the text NASA ADS  
• Semenov, D., Henning, Th., Helling, Ch., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 In the text NASA ADS  
• Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 In the text NASA ADS  
• Snyder, L., Veal, J., Woodney, L., et al. 2001, ApJ, 121, 1147 In the text  
• Spalding, D. B. 1971, Chem. Eng. Sci., 26, 95 In the text  
• Strehmel, K. & Weiner, R. 1995, Numerik gewöhnlicher Differentialgleichungen (Stuttgart: Teubner) In the text  
• Toomre, A. 1964, ApJ, 139, 1217 In the text NASA ADS  
• Taylor, S., Morata, O., & Williams, D. 1996, A&A, 313, 269 In the text NASA ADS  
• Umebayashi, T., & Nakano, T. 1980, Publ. Astron. Soc. Jpn, 32, 405 In the text NASA ADS  
• Umebayashi, T., & Nakano, T. 1981, Publ. Astron. Soc. Jpn, 33, 617 In the text NASA ADS  
• Unsöld, A. 1968, Physik der Sternatmosphären (Berlin: Springer) In the text  
• Wehrstedt, M., & Gail, H.-P. 2002, A&A, 385, 181 NASA ADS  
• Willacy, K., Klahr, H., Millar, T., & Henning, Th. 1998, A&A, 388, 995 In the text  
• Willacy, K., & Langer, W. 2000, ApJ, 544, 903 In the text NASA ADS  
• Wutz, M., Adam, H., & Walcher, W. 1992, Theorie und Praxis der Vakuumtechnik, Vieweg, Braunschweig In the text  
• Xie, T., Allen, M., & Langer, W. 1995, ApJ, 440, 674 In the text NASA ADS