<figure><img src="image_1.png" /></figure>

HAL
open science

# The vertical structure of T Tauri accretion disks. II.
## Physical conditions in the disk

Fabien Malbet, Régis Lachaume, Jean-Louis Monin

► To cite this version:

Fabien Malbet, Régis Lachaume, Jean-Louis Monin. The vertical structure of T Tauri accretion disks. II. Physical conditions in the disk. Astronomy & Astrophysics - A&A, 2001, 379 (2), pp.515-528. 10.1051/0004-6361:20011339 . hal-05333779

HAL Id: hal-05333779

https://hal.science/hal-05333779v1

Submitted on 28 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The vertical structure of T Tauri accretion disks

## II. Physical conditions in the disk

F. Malbet¹, R. Lachaume¹, and J.-L. Monin¹,²

¹ Laboratoire d'Astrophysique UMR UJF-CNRS 5571, Observatoire de Grenoble, Université Joseph Fourier, BP 53, 38041 Grenoble Cedex 9, France

² Institut Universitaire de France

Received 18 May 2001 / Accepted 11 September 2001

**Abstract.** We present a self-consistent analytical model for the computation of the physical conditions in a steady quasi-Keplerian accretion disk. The method, based on the thin disk approximation, considers the disk as concentric cylinders in which we treat the vertical transfer as in a plane-parallel medium. The formalism generalizes a work by Hubeny (1990), linking the disk temperature distribution to the local energy dissipation and leads to analytical formulae for the temperature distribution which help to understand the behaviour of the radiation propagated inside the disks. One of the main features of our new model is that it can take into account many heating sources. We apply the method first to two sources: viscous dissipation and stellar irradiation. We show that other heating sources like horizontal transfer or irradiation from the ambient medium can also be taken into account. Using the analytical formulation in the case of a modified Shakura & Sunyaev radial distribution that allow the accretion rate to be partly self-similar in the inner region, and, for an $\alpha$ and $\beta$ prescription of the viscosity, we obtain two-dimensional maps of the temperature, pressure and density in the close environment of low mass young stars. We use these maps to derive the observational properties of the disks such as spectral energy distributions, high resolution spatial images or visibilities in order to underline their different behaviours under different input models.

**Key words.** accretion, accretion disk – radiative transfer – stars: circumstellar matter – stars: pre-main sequence

## 1. Introduction

Since the initial models of viscous accretion disks of Shakura & Sunyaev (1973) and Lynden-Bell & Pringle (1974), many efforts have focused on the description of the close environment of T Tauri stars (hereafter TTS). In order to reproduce high angular resolution images that became available soon after 1990, various authors have published models of scattered light in large circumstellar structures (Lazareff et al. 1990; Whitney & Hartmann 1992). In the last five years, more precise images of edge-on disks have been obtained by the Hubble Space Telescope (HST) and infrared adaptive optics (Burrows et al. 1996; Stapelfeldt et al. 1998; Monin & Bouvier 2000), allowing us to refine the models and the physical parameters used to describe the disks. Most of the time, when interpreting these images, the authors use ad-hoc power laws to extract the physical conditions in the disks, like surface density or the disk scale height. Indeed, these images show that the

large-scale flaring disks have a complex vertical structure that must be taken into account in the interpretation of their flat SED.

Up to now, we have access to only a very few constraints on the very central parts of circumstellar disks, since current instruments do not allow us to directly observe within 10 AU of the central object. However, the physics of the inner part of the disks is not expected to be completely different from that in their outer parts, but our knowledge of these central regions of circumstellar disks is restricted to scarce interferometric data and indirect measurements of the inner phenomena via the spectral energy distribution (hereafter SED), the high energy UV and X-ray emission of disks or the measurement of magnetic fields.

The first models used quasi-Keplerian steady accretion disks, expected to be geometrically thin for reasonable values of the accretion rate. They describe the disk as infinitely flat, predicting that the emergent flux has a spectral energy distribution proportional to the power $-4/3$ of the wavelength in the infrared. However,

<footnote>Send offprint requests to: F. Malbet,
e-mail: Fabien.Malbet@obs.ujf-grenoble.fr</footnote>
Rydgren & Zak (1987) showed that most TTS present flatter spectral energy distributions, i.e. decreasing less steeply or even remaining constant over a large domain in the infrared. This discrepancy led Adams et al. (1988) to assume that the radial temperature distribution in T Tauri disks is also flatter than in standard accretion disks: if one assumes that the disk thickness flares then the disk intercepts more light from the star at large distances and leads to a variety of SED slopes.

Modeling of the vertical structure of the disks and the associated radiative transfer has been developed by Bell & Lin (1994) focusing their attention on the vertical structure of disks in FU Orionis systems where they explained a possible origin of outbursts. Chiang & Goldreich (1997, 1999) used a semi-analytical model of the circumstellar structure to explain the SEDs of different types of young stars. Recently, D’Alessio et al. (1998, 1999) numerically resolved a complete set of equations to compute the spatial distribution of the temperature and the density that can be constrained by the observed SEDs and two-dimensional images.

Up to now, the models have included various approximations: (i) the magnetic field is ignored; (ii) the viscosity is assumed to depend only on local conditions and (iii) the viscosity is expressed using a simple ad-hoc law. On the other hand, when authors take into account the magnetic field (Shu et al. 1994; Ferreira & Pelletier 1995; Casse & Ferreira 2000a; Casse & Ferreira 2000b), they assume that these magnetic fields do not affect the overall disk structure. As another example, Terquem et al. (1998) and Terquem (1998) showed that gravitational instabilities in the disks may produce waves dissipating viscous energy far from where they are created. Integrating all these effects in a unique model appears to be a difficult task.

Our goal is to develop a radiative model simple enough that it could later on include as many detailed physical processes of the energy transport in the disks as possible, and/or the interaction between the circumstellar matter and the magnetic field. In this paper, we present a new method to process the radiative transfer in circumstellar disks, following the pioneering work by Hubeny (1990) and an initial application in the case of TTS disks by Malbet (1991, 1992). Our results are not fundamentally different from previous models but thanks to the analytical expression of the vertical temperature distribution that we develop, we are able to identify the actual origin of the temperature rise in the upper part of the disk photosphere (1991). We propose that this formalism can be the starting point for more elaborate models, including magnetic field or viscous transport. As a first step, our new self-consistent analytical model can be used to constrain the physical environment of young stars using images from HST or ground-based large telescopes, as well as millimeter radio-interferometer maps and forthcoming data from large optical interferometers, and especially the ESO VLTI project (Malbet et al. 2000).

Section 2 presents the analytical model derived from hydrostatic and radiative equilibrium, as first described

by Hubeny (1990), and extended here to the illumination by an external source. Section 3 describes the application of this model to the specific case of the circumstellar environment of young stellar objects, especially in the case where the accretion rate is not constant at all radii. To illustrate the model and its application, we present and discuss in Sect. 4 some results in the case where the heating by the central source can be neglected and we focus our attention on the production of observable quantities like SEDs, images, and interferometric data including visibility amplitudes and closure phases.

## 2. Analytical model

In this section, we detail the analytical model for the calculation of the spatial distribution of quantities like the temperature or the density. We use the notations of Hubeny (1990) as much as possible, even if we present the formalism in a slightly different way. The first 3 sections (2.1, 2.2 and 2.3) are short reminders of the main equations of the problem. The last one (2.4) is the generalization of this formalism to external sources of radiation.

The principle for computing the vertical structure of T Tauri disks is to consider the disk as a set of concentric cylinders of infinitesimal width. Each cylinder is itself regarded as a plane-parallel atmosphere of mixed dust and gas layers. The basic assumption is that the radial radiative transfer of energy is smaller than the vertical transfer, which is true only in the thin disk approximation (Pringle 1981). However, we can still take into account the horizontal transfer as a perturbation for disks of moderate flaring (see in Sect. 2.4.4). We solve the coupled equations of hydrostatic equilibrium and the radiative transfer in each cylinder in order to compute the temperature and density distributions along the z-axis at a given radius $r$. In the following sections, when we focus our attention on the vertical structure of the cylinder of radius $r$, we drop the $r$ coordinate unless otherwise specified. We also use the mass column coordinate $m$ in place of the height coordinate $z$:

$$m(z) = \int_{z}^{\infty} \rho(z) dz, \qquad (1)$$

where $M = m(0) = \Sigma/2$ is half the surface density of the disk at radius $r$.

### 2.1. Initialization stage from the radial structure

For the initial step, we assume a standard geometrically thin disk with a uniform vertical structure at each radius. To fully solve the radial structure, we need:

- the mass per area unit $\Sigma(r)$

- the energy dissipation per area unit $F(r)$

- the opacity law $\chi(T, \rho)$ of the disk material

- the equation of state of the material, $\rho = S(P, T)$

- the local vertical acceleration field $G_z(r, z)$, generally inferred from the gravity
where the temperature, density, scale height, pressure and optical thickness laws, respectively $T_D$, $P_D$, $\rho_D$, $h_D$ and $\tau_D$, are given by the coupled equations:

$$ 
\begin{aligned}
T_D^4(r) &= \frac{F(r)/2}{1 - e^{-\tau_D(r)}} \\
P_D(r) &= -G_z(r, h_D(r)) \Sigma(r) \\
\rho_D(r) &= S(P_D(r), T_D(r)) \\
h_D(r) &= \Sigma(r)/\rho_D(r) \\
\tau_D(r) &= \Sigma(r) \chi(T_D(r), \rho_D(r)).
\end{aligned} 
\quad (2) $$

In the optically thin case, an iterative computation is required. One should notice that the mass, energy dissipation and acceleration field distributions do not need to be independent of the temperature or density laws; in this case $\Sigma(r)$, $F(r)$ and $G_z(r, z)$ are refined at each iteration.

A first refinement is to replace the uniform distribution $\rho_D(r)$ by the isothermal density distribution. The resulting radial structure is then the structure given by the so-called standard model (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974).

Another improvement consists of replacing the effective temperature by the temperature of the equatorial layer located at $\tau_D/2$. As in a stellar atmosphere, we substitute for the first equation of the set of Eqs. (2) the following equation:

$$ T_D^4(r) = \frac{3/4 F(r) (\tau_D(r) + 2/3)}{1 - e^{-\tau_D(r)}}. \quad (3) $$

## 2.2. Hydrostatic equilibrium – Density distribution

For a given vertical temperature law, we compute the density distribution from the pressure given by the hydrostatic equilibrium. As the disk is dominated by gravitation of the central star, the partial derivative of $G_z$ is a constant in the approximation $z \ll r$, and we use a differentiated hydrostatic equilibrium equation for convenience:

$$ \frac{d^2 P}{dm^2} = - \frac{1}{S(P, T)} \frac{\partial G_z}{\partial z}. \quad (4) $$

The second order equation Eq. (4) is solved with two following boundary conditions:

1. in the disk mid-plane: the first derivative of $P$ is zero by symmetry;

2. in the outer boundary: $P$ tends to a constant low value, the pressure in the interstellar medium (hereafter ISM). The choice of this limit has almost no influence on the disk structure.

## 2.3. Radiative equilibrium – Temperature distribution

To complete the calculations of the distributions of the physical quantities, we have to connect the local temperature with the heating sources, remote or local. This is the goal of the radiative transfer described in this section.

Theoretically speaking, the radiative transfer equations can be applied to any kind of radiation. However,

common and well-known approximations make them easier to solve provided the radiation is thermalized or isotropic. For instance we can approximate the flux-weighted mean opacity by the Rosseland opacity if the radiation is thermalized with the medium. The Eddington closure relations apply well for isotropic radiation.

We therefore treat separately radiation emitted (from internal heating sources or from external light absorbed by the disk) which is mainly isotropic and thermalized, and the radiation coming from external sources. The latter can also be split into two terms, the attenuated radiation coming from its external source in a certain known direction, and the radiation that has been scattered one or several times which should be mainly isotropic. In the first part of this section, following Hubeny we neglect the radiation by external sources. Later in Sect. 2.4, we will take into account correctly the influence of external radiations.

### 2.3.1. Energy dissipation in each layer

We first express the temperature as a function of the local production of energy $u$, such as viscous heating, stellar irradiation, disk backwarming, heating by a UV or X radiation field, or the contribution of the energy horizontal transfer. Following the usual notations (Mihalas 1978), we link $u$ to the source function $S$ and the diffuse radiation intensity $I$:

$$ 4\pi \int_{0}^{\infty} \chi_{\nu}(m) [S_{\nu}(m) - J_{\nu}(m)] \, d\nu = \rho(m)u(m) \quad (5) $$

$$ \text{with } J_{\nu}(m) = \int_{-1}^{1} I_{\nu}(m, \mu) \, d\mu. \quad (6) $$

The conservation of the radiative energy means that in an atmosphere layer, the radiated energy, i.e. the emission of the layer, minus the radiation coming from the neighbouring layers $J_{\nu}(z)$, is exactly equal to the dissipated energy in the layer:

$$ S_{\nu} = \frac{\kappa_{\nu}}{\chi_{\nu}} B_{\nu} + \left(1 - \frac{\kappa_{\nu}}{\chi_{\nu}}\right) J_{\nu}. \quad (7) $$

Integrating over frequencies, one derives from Eq. (5):

$$ \kappa_B B = \kappa_J J + \frac{u}{4\pi} \quad (8) $$

where $\kappa_B$ and $\kappa_J$ are respectively the Planck absorption mean and the usual absorption mean, defined by:

$$ \kappa_B = \int_0^\infty (\kappa_\nu / \rho) B_\nu \, d\nu / B \quad (9) $$

$$ \kappa_J = \int_0^\infty (\kappa_\nu / \rho) J_\nu \, d\nu / J. \quad (10) $$

Since $\pi B = \sigma_B T^4$, the vertical temperature distribution is, by using Eq. (8):

$$ T^4 = \frac{\pi \kappa_J}{\sigma_B \kappa_B} \left[ J + \frac{u}{4\pi \kappa_J} \right]. \quad (11) $$
In other words, the temperature in a layer is the sum of two contributions: (i) a term proportional to the radiative intensity $J$ transferred by adjacent layers and (ii) a local dissipation term proportional to $u$.

### 2.3.2. Radiative transfer of the diffuse intensity

We solve the transfer of the diffuse intensity $I$ to determine $J$. The equation of transfer and its two first moments are:

$$ \mu \frac{\partial I_{\nu}(z, \mu)}{\partial z} = \chi_{\nu}(z) [S_{\nu}(z) - I_{\nu}(z, \mu)], \quad (12) $$

$$ \frac{\mathrm{d} H_{\nu}(z)}{\mathrm{d} z} = \chi_{\nu}(z) [S_{\nu}(z) - J_{\nu}(z)], \quad (13) $$

$$ \frac{\mathrm{d} K_{\nu}(z)}{\mathrm{d} z} = -\chi_{\nu}(z) H_{\nu}(z) \quad (14) $$

where $J_ν$, $H_ν$ and $K_ν$ are respectively the zeroth, first and second moments of $I_ν$.

By integrating over frequencies one gets:

$$ \frac{\mathrm{d} H}{\mathrm{d} m} = -\frac{u}{4\pi}, \quad (15) $$

$$ \frac{\mathrm{d} K}{\mathrm{d} m} = \chi_H H. \quad (16) $$

In this equation, the flux-averaged opacity $\chi_H$ is defined by:

$$ \chi_H = \int_0^\infty (\chi_\nu / \rho) H_\nu \, d\nu / H. \quad (17) $$

When integrating the first moment equation Eq. (15) over $m$, one gets:

$$ H(m) = H(0) - U M \theta(m) / 4\pi \quad (18) $$

$$ \text{with } U M = \int_{0}^{M} u(\zeta) \, d\zeta \quad (19) $$

$$ \text{and } \theta(m) = \frac{1}{UM} \int_{0}^{m} u(\zeta) \, d\zeta, \quad (20) $$

where $U$ is the total energy dissipation in the upper atmosphere of the disk per surface unit, and $\theta$ the distribution function of the energy dissipation in this atmosphere.

The symmetry conditions in the equatorial plane for the flux $H(m)$ implies that $H(M) = 0$, that is $\theta(M) = 1$. We therefore get a relation between $H(0)$ and the energy flux dissipated in all atmospheric layers $U$:

$$ H(0) = U M / 4\pi. \quad (21) $$

Unlike in a stellar atmosphere, energy is locally produced so that the flux is not constant.

Integrating the second moment equation Eq. (16) over $m$, one derives:

$$ K = K(0) + H(0) (\tau_H - \Delta\tau_H) \quad (22) $$

$$ \text{with } \tau_H(m) = \int_0^m \chi_H(\zeta) \, d\zeta \quad (23) $$

$$ \text{and } \Delta\tau_H(m) = \int_0^m \chi_H(\zeta)\theta(\zeta) \, d\zeta \quad (24) $$

where $\tau_H$ and $\Delta\tau_H$ are respectively the flux weighted and the “energy dissipation weighted” mean optical depths. In a stellar atmosphere, one has $K = K(0) + H(0)\tau_H$. Here, $\Delta\tau_H$ represents a correction to the flux-weighted mean optical depth due to the decrease of the flux toward the mid-plane of the disk. In the large depth approximation, we have: $\mathrm{d}T^4/\mathrm{d}\tau_H \propto H$. Since in a stellar atmosphere, the flux $H$ is constant, $T^4 \propto H\tau_H$, whereas in a disk $H$ decreases with $\tau_H$ until $T^4 \le H(0)\tau_H$, hence the corrective term $\Delta\tau_H$.

We then introduce the Eddington factors,

$$ f_K = \frac{K}{J} \quad \text{and} \quad f_H = \frac{H(0)}{J(0)}. \quad (25) $$

In the large depth approximation, met in thick parts of the disk, they are respectively equal to 1/3 and 1/2.

### 2.3.3. Formal solution

Combining Eqs. (11), (22) and (25), we get the temperature distribution:

$$ T^4 = \frac{\kappa_J T_{\text{eff}}^4}{4\kappa_B f_K} \left[ \left( \tau_H - \Delta\tau_H + \frac{f_K(0)}{f_H} \right) + \frac{f_K}{M\kappa_J} \frac{u}{U} \right] \quad (26) $$

with $\sigma_B T_{\text{eff}}^4 = UM$.

In the case of strict thermal equilibrium, the intensity-weighted and flux-weighted mean opacities $\kappa_J$ and $\kappa_H$ can be substituted by the Planck and Rosseland mean opacities $\kappa_B$ and $\kappa_R$.

Equation (26) means that the vertical distribution of the temperature in a disk is similar to its distribution in a stellar atmosphere:

$$ T^4(\tau) = \frac{\kappa_J(\tau)}{4\kappa_B(\tau)f_K(\tau)} T_{\text{eff}}^4 \left[ \tau + \frac{f_K(0)}{f_H} \right], \quad (27) $$

with two main differences:

(i) A local term, proportional to the energy dissipation $u$ is added. In a stellar atmosphere no energy dissipation occurs in the intermediate layers. The factor proportional to the density of energy per unit mass $u$, in Eq. (26) is also inversely proportional to the optical depth of the disk, $\kappa_J M$. This factor is therefore dominant in zones where the disk is thin;

(ii) the optical depth $\tau$ becomes $\tau_H - \Delta\tau_H$. The term $\Delta\tau_H$ is positive but cannot exceed $\tau_H$ depending where the energy dissipation occurs in the vertical structure. The two extreme situations are when the dissipation function is a Dirac function located either in the central plane of the disk, or at the top of the photosphere. In the first situation, $\Delta\tau_H = 0$ everywhere, whereas in the second one $\Delta\tau_H = \tau_H$. Above the layers where most dissipation occurs, the correction $\Delta\tau_H$ can be ignored and the disk behaves like a stellar atmosphere: there is no local energy source and all radiation comes from deeper layers. In the central layers, $\theta(m) \approx 1$ and $\tau_H(m) - \Delta\tau_H(m)$ is constant: most energy dissipation occurs in upper layers, and, in other terms
there is almost no radiation coming from deeper layers. Therefore the temperature is more or less constant. In the particular case where the energy is dissipated only in the central plane $\tau_H - \Delta\tau_H = \tau_H$ and the disk photosphere behaves like a stellar atmosphere.

### 2.3.4. Several sources of energy dissipation

In the case where there are several sources of energy dissipation releasing the energy $u_k(m)$ in the layer $m$, the formal solution given in Eq. (26) can be split as the sum of the temperature distributions (see Appendix A):

$$T^4(m) = \sum_k T_k^4(m), \qquad (28)$$

with the temperature distribution of source $k$,

$$T_k^4(m) = \frac{\kappa_J(T_{\text{eff}})_k^4}{4\kappa_B f_K} \times \left[ \left( \tau_H - (\Delta\tau_H)_k + \frac{f_K(0)}{f_H} \right) + \frac{f_K}{M\kappa_J} \frac{u_k}{U_k} \right], \quad (29)$$

where the effective temperature of the $k$th source is defined by:

$$\sigma_B(T_{\text{eff}})_k^4 = \int_0^M u_k(\zeta) \, d\zeta = U_k M. \qquad (30)$$

The quantity $(\Delta\tau_H)_k$ is defined by Eq. (24) where the vertical distribution function $\theta(m)$ is replaced by the corresponding $k$th energy source,

$$\theta_k(m) = \frac{1}{U_k M} \int_0^m u_k(\zeta) \, d\zeta. \qquad (31)$$

We see that when dealing with several sources of energy dissipation, one can treat them separately from the analytical point of view (see Sect. 2.4). However in the full modelization, the field of radiation which defines $J(m)$, $H(m)$, $f_K(m)$ and $f_H$ must be computed globally.

### 2.4. Influence of external radiative sources

Hubeny (1990) already took into account the effect of the irradiation by the central star (see also Hubeny 1991). However, when addressing the influence of the radiation from external sources, we can use Eqs. (28) and (29) taking into account as an additional source of energy dissipation the reprocessing of the incoming radiation. To compute the reprocessed energy $u_A$ and its influence on $(\Delta\tau_H)_A$, we need to solve the radiative transfer of the direct and scattered incoming radiation.

The flux coming from a central source is highly directional and cannot be processed as isotropic radiation. For convenience, one may split the radiation field into several components, and consider separate partial transfer equations for the individual components. Some of the radiation components may be considered as a source of additional

heating. Only the component containing the thermal emission term, $\kappa_\nu B_\nu$ can be treated with the formalism developed in the previous section. The mean opacities $\kappa_J$, $\kappa_B$, $\chi_H$, and the mean Eddington factors $f_K$ and $f_H$ are then defined as appropriate averages over the moments of this particular component of the radiation field.

Here, we process separately (i) the radiation from internal heating sources or from external light absorbed and reprocessed in the disk, a radiation which is mainly isotropic and thermalized, and (ii) the radiation coming from external sources which is not absorbed. In this latter case, we again split the incoming radiation into two terms, one attenuated coming directly from its external source in a given known direction, and one that has been scattered once or several times and that should be mainly isotropic.

We use the following notations:

1. $I_\nu$: radiation emitted by the disk itself. (All thermal flux related quantities are written without superscript);

2. $I_\nu^s$: radiation coming from external sources and scattered in the disk;

3. $I_\nu^0$: attenuated radiation coming from external sources

4. $u_\nu$: the energy dissipated locally (e.g. viscosity process);

5. $u_A$: the energy coming from the incident photons absorbed in the disk and reprocessed.

This approach is quite similar to Chandrasekhar's one (Chandrasekhar 1960) and was used by (1991). It consists of writing the equations of transfer for the total intensity $I_\nu$ and then deriving equations involving non-local terms $I_\nu^0$ and $I_\nu^s$. Since the radiative transfer of $I_\nu$ and non-local terms are treated separately, the coupling between them can be seen either as a loss of energy for the non-local radiation, or as a gain of energy for $I_\nu$. The stellar flux can now be seen as energy dissipation similar to that due to viscosity.

#### 2.4.1. Energy reprocessed by the disk

In the presence of an external source of radiation, the energy dissipation term $u_A$ linked to the reprocessing of the radiation coming from external sources and absorbed in the disk is directly linked to the specific intensity of the incoming radiation by:

$$\rho(m)u_A(m) = 4\pi \int_{-1}^{+1} \int_0^{+\infty} \kappa_\nu (I_\nu^0(m, \mu) + I_\nu^s(m, \mu)) \, d\nu \, d\mu. \qquad (32)$$

With $J^0$ and $J^s$ the frequency-integrated zeroth order moments of $I_\nu^0$ and $I_\nu^s$,

$$J^0 = \int_0^{+\infty} J_\nu^0 \, d\nu \quad \text{with} \quad J_\nu^0 = \int_{-1}^{+1} I_\nu^0(\mu) \, d\mu, \qquad (33)$$

$$J^s = \int_0^{+\infty} J_\nu^s \, d\nu \quad \text{with} \quad J_\nu^s = \int_{-1}^{+1} I_\nu^s(\mu) \, d\mu, \qquad (34)$$
and, the intensity weighted opacities,

$$ \kappa_J^0 = \int_0^{+\infty} (\kappa_\nu/\rho) J_\nu^0 d\nu / J^0, \quad (35) $$

$$ \kappa_J^s = \int_0^{+\infty} (\kappa_\nu/\rho) J_\nu^s d\nu / J^s. \quad (36) $$

u_A can be expressed as:

$$ u_A = 4\pi(\kappa_J^0 J^0 + \kappa_J^s J^s). \quad (37) $$

The energy absorbed u_A is therefore proportional to the mean incident flux and to the absorption coefficient of the medium.

### 2.4.2. Radiative transfer of the incoming radiation

We study the transfer of the radiative intensities $I^0$ and $I^s$ upon which our knowledge of u_A and $(\Delta\tau_H)_A$ depends. The equations of transfer for the incoming radiation within the frame of coherent isotropic scattering are:

$$ \mu \frac{\partial I_\nu^0}{\partial z}(z, \mu) = -\chi_\nu I_\nu^0(z, \mu), \quad (38) $$

$$ \mu \frac{\partial I_\nu^s}{\partial z}(z, \mu) = -\chi_\nu I_\nu^s(z, \mu) + \sigma_\nu J_\nu^s(z, \mu) + \sigma_\nu J_\nu^0(z, \mu), \quad (39) $$

where $\sigma_\nu J_\nu^0(z, \mu)$ accounts for the scattering of the incoming radiation and $\sigma_\nu J_\nu^s(z, \mu)$ for multiple scattering.

We define the direction and frequency averaged opacities,

$$ \chi_J^0 = \int_0^{+\infty} (\chi_\nu/\rho) J_\nu^0 d\nu / J^0; \quad (40) $$

$$ \chi_H^0 = \int_0^{+\infty} (\chi_\nu/\rho) H_\nu^0 d\nu / H^0; \quad (41) $$

$$ \sigma_J^0 = \int_0^{+\infty} (\sigma_\nu/\rho) J_\nu^0 d\nu / J^0; \quad (42) $$

$$ \chi_H^s = \int_0^{+\infty} (\chi_\nu/\rho) H_\nu^s d\nu / H^s; \quad (43) $$

where $J_\nu^0$, $J_\nu^s$, and, $H_\nu^0$, $H_\nu^s$ are respectively the zeroth and first moments of $I_\nu^0$ and $I_\nu^s$. Therefore the frequency-integrated zeroth, first and second moments of $I_\nu^0$ and $I_\nu^s$, namely $J^0$, $H^0$, $K^0$ for direct incoming radiation and $J^s$, $H^s$, $K^s$ for scattered incoming radiation, are linked by the following equations:

$$ \frac{d H^0}{dm} = -\chi_J^0 J^0 \quad \text{and} \quad \frac{d K^0}{dm} = -\chi_H^0 H^0 \quad (44) $$

$$ \frac{d H^s}{dm} = -\kappa_J^s J^s + \sigma_J^0 J^0 \quad \text{and} \quad \frac{d K^s}{dm} = -\chi_H^s H^s. \quad (45) $$

Therefore, the reprocessed energy,

$$ u_A(m) = -4\pi \left( \frac{d H^0}{dm} + \frac{d H^s}{dm} \right), \quad (46) $$

leads to

$$ U_A M = 4\pi (H^0(0) + H^s(0)) \quad (47) $$

$$ \theta_A(m) = 1 - \frac{H^0(m) + H^s(m)}{H^0(0) + H^s(0)} \quad (48) $$

assuming $H^s(M) = 0$ and $H^0(M) = 0$, i.e. assuming that the equatorial plane is a plane of symmetry for the source of external radiation. This is truly the case for a central star or a uniform ambiant medium.

Using the second moment of the incoming and reprocessed radiative equations, we derive $(\Delta\tau_H)_A$,

$$ (\Delta\tau_H)_A = \tau_H + \frac{\omega^0 [K^0 - K^0(0)] + \omega^s [K^s - K^s(0)]}{H^0(0) + H^s(0)} \quad (49) $$

where we define respectively $\omega^0$ and $\omega^s$, the mean ratios of $\chi_H$ to $\chi_H^0$ and $\chi_H^s$:

$$ \int_0^m \chi_H(\zeta) H^0(\zeta) d\zeta = \omega^0(m) \int_0^m \chi_H^0(\zeta) H^0(\zeta) d\zeta \quad (50) $$

$$ \int_0^m \chi_H(\zeta) H^s(\zeta) d\zeta = \omega^s(m) \int_0^m \chi_H^s(\zeta) H^s(\zeta) d\zeta. \quad (51) $$

If the radiation is absorbed in a geometrically thin layer, where the temperature and density are more or less constant, then $\omega^0 \approx \chi_H/\chi_H^0$. It is the ratio of the absorption of reprocessed radiation (generally in the thermal infrared) to the absorption of incoming radiation (often in the visible, or UV, or X-ray domain). In the case of a hot source like stellar irradiation or a UV field, this ratio is expected to be less than unity.

Finally the contribution of the incoming radiation source to the temperature distribution is:

$$ T_A^4(m) = \frac{\kappa_J}{4\kappa_B f_K} (T_{\text{eff}})_A^4 \times \left[ \frac{\omega^0 [K^0(0) - K^0] + \omega^s [K^s(0) - K^s]}{H^0(0) + H^s(0)} + \frac{f_K(0)}{f_H} + \frac{f_K}{M\kappa_J} \frac{\kappa_J^0 J^0 + \kappa_J^s J^s}{H^0(0) + H^s(0)} \right] \quad (52) $$

We show in Appendix B how to compute $J^0(m)$, $J^s(m)$, $H^0(m)$, and $H^s(m)$ in the case of an external point-like source where the irradiation from the central source is highly collimated. As demonstrated in Sect. 2.3.4, the temperature vertical structure of a sum of point-like sources is the sum of the temperatures computed for each source to the fourth power.

### 2.4.3. Interpretation

In this section, we discuss in detail the physical interpretation of the complex Eq. (52). Since we know how to add the contributions of various energy sources by using Eqs. (28) and (29), we can restrict ourselves to a point-like source located at infinity without losing generality. An extended source will then be processed as the sum of point-like sources.

The dependence of the central temperature to the incoming stellar radiation is relatively simple. Since the topmost layer emits half upward to the outside and half downward toward the inner layers, the latter ones receive half of the emission. Their effective temperature is then $T_A/2^{1/4}$. Since there is no flux inside (for $\tau \gg 1$), the temperature
is constant and equal to this temperature. This intuitive result is compatible with Eq. (52). As a matter of fact, if $\mu^0$ is the cosine angle of the incoming radiation and since $K = \mu^0 H$, the first term in the square brackets is of the order of $\omega^0 \mu^0$. For a slanted incident stellar radiation, we have $\mu^0 \ll 1$ and therefore $\omega^0 \mu^0 \ll 1$. Deep inside the disk, the third term, evanescent since the incoming radiation is absorbed, is also negligible. The brackets reduce to $f_K(0)/f_H \approx 2/3$. The central temperature is therefore smaller than the effective temperature $T_A$ by a factor $\approx 2^{-1/4}$. In conclusion, half of the absorbed energy is irradiated by the hot topmost layer and another half by the cooler inner layers.

At the surface, the third term in the square brackets reduces to $1/(\mu^0 \omega^0)$ by using $J^0 = H/\mu^0$ and is dominant. The temperature is therefore much higher than the effective temperature. This fact is explained in Malbet & Bertout (1991): the topmost layer, vertically optically thin but radially optically thick, absorbs almost all stellar light. This layer is superheated by a factor of $(1/\mu^0)^{1/4}$. Chiang & Goldreich (1997) refined this view by multiplying the term $1/\mu^0$ by $1/\omega^0$. The first factor is a pure geometrical effect of the slanted incidence whereas the second one comes from the difference between the opacities for the visible incoming radiation and the infrared reprocessed radiation.

### 2.4.4. Overview of other heating sources

The flux radiation from the ISM can be processed on the same basis as the stellar irradiation. The main difference is that the ISM does not irradiate in a particular direction. If we assume that it consists of an isotropic black body radiation of characteristic temperature $T_{ISM}$, we can apply Eq. (61) within the two-stream approximation (i.e. with $\mu^* = 1/\sqrt{3}$ using the definition given in Sect. 3.2.2) and we can demonstrate that the contribution of the ISM radiation is a constant. In fact, if the ISM is the unique source of irradiation, the disk would be isothermal.

All other radiative sources can be treated exactly as the central star. If they are extended, the previous results derived for point-like sources are integrated over angles.

The radial extent of the disk is much larger than the vertical one. We can therefore apply the approximations of large depth to the horizontal transfer. This hypothesis, combined with the LTE approximation, allows us to write the horizontal flux $H_r$ as:

$$H_r = -\frac{4}{3\pi(\rho\kappa_R)}T^3\frac{\partial T}{\partial r} \qquad (53)$$

The local energy gain for the material between $r$ and $r+dr$ is $4\pi(H_r(r+dr)-H_r(r))$. One then derives the heating per unit of mass of material that allow us to self-consistently compute the correction to the vertical transfer, in the case $H_r \ll H$:

$$u_H = \frac{4\pi}{\rho} \frac{\partial H_r}{\partial r} \qquad (54)$$

## 3. Application to T Tauri disks

The two main heating sources in the case of T Tauri disks are the viscous dissipation and the stellar light absorption. We apply our general equations in this frame.

### 3.1. Approximations

All the terms $\tau_H$, $\kappa_H$, etc., are dependent on the structure. Moreover $f_K$ and $f_H$, are not known a priori. Therefore we will seek a self-consistent iterative method able to determine these quantities.

We now make two approximations concerning the opacities: the diffuse intensity is a Planck distribution at the local temperature and the radiation coming from a external source is also a Planck distribution at the temperature of the source $T^*$. We then derive:

- $\kappa_J = \kappa_B$, the Planck mean opacity.

- $\kappa_J^0 = \kappa_J^s = \kappa_B^*$, the mean opacity for a black body radiation at $T^*$ through a medium at temperature $T$.

- $\kappa_H = \kappa_R$, the Rosseland mean opacity.

Since the scope of this paper is to present analytical solutions and simple numerical applications, we use the gray Rosseland opacities given by Bell & Lin (1994). We also use the gray approximation $\kappa_B = \kappa_B^* = \kappa_R$ and ignore diffusion. Therefore $\omega^s = \omega^0$ that we now call $\omega$.

The $f_H$ and $f_K$ values will be close to 1/2 and 1/3. These approximations are only valid at large depths, but we refine these values in optically thin zones at each iteration, as well as for $\omega$.

### 3.2. Vertical temperature distribution

#### 3.2.1. Viscous heating at each layer in the disk

We consider here an energy density per mass unit which depends on local conditions. The energy dissipation per mass unit due to viscosity is given by Franck et al. (1985):

$$u_V(m) = \nu_{\text{visc}}(m) \left( r \frac{d\Omega}{dr} \right)^2 \qquad (55)$$

In the case of a thin Keplerian disk we get:

$$u_V(m) = \frac{9}{4} \Omega^2 \nu_{\text{visc}}(m) \qquad (56)$$

Since $u_V(m)$ depends vertically only on $\nu_{\text{visc}}(m)$, the $\theta_V$ function depends only on the vertical structure of the viscosity $\nu_{\text{visc}}$ and of its mass-average $\bar{\nu}_{\text{visc}}$.

Up to now there has been no exact description of the viscosity induced by the turbulence. In this paper, we have considered two different types of viscosity:

1. The Shakura-Sunyaev prescription stating that the viscosity is proportional to the local sound speed and local height scale:

$$\nu_{\text{visc}} = \alpha c_s h \propto T. \qquad (57)$$
The $\theta_V$ function is therefore the distribution function of the temperature. The term $\tau_H - (\Delta\tau_H)_V$, by definition smaller than $\tau_H$, can be much smaller if the temperature peaks at the surface of the photosphere;

2. The $\beta$-prescription, derived from laboratory experiments, and applied to accretion disks by Huré et al. (2001), states that

$$ \nu_{\text{visc}} = \beta \Omega r^2 \quad (58) $$

where $\beta$ is a constant ranging from $10^{-5}$ to $10^{-3}$. Here the viscosity is uniform along the vertical axis and does not depend on the temperature distribution of the disk.

### 3.2.2. Heating by the central stellar source

With the approximations of Sect. 3.1, there is no reason to distinguish the two specific intensities $I^0$ and $I^s$. If we call $I^* = I^0 + I^s$, then we find an analytical solution very similar to Eq. (28) of from Malbet & Bertout (1991):

$$ T_A^4 = \frac{\pi \kappa_J H^*(0)}{f_K \kappa_B \sigma_B} \left[ \omega^* \frac{K^*(0) - K^*}{H^*(0)} + \frac{f_K(0)}{f_H} + \frac{f_K \kappa_J^* J^*}{\kappa_J H^*(0)} \right] \quad (59) $$

where the quantities linked to the star incoming radiation are superscripted with *. The only difference is the presence of $\omega^*$ which is not a priori equal to 1 as noted by D’Alessio et al. (1998).

We derive $J^*(0)$ and $H^*(0)$ from the geometrical properties of the system using the approximation described by Ruden & Pollack (1991). The star is seen by the disk as a point-like source. It is valid beyond a few stellar radii; in the inner parts this is not the case, nevertheless re-processing is seldom dominant in this region (see Fig. 3 of D’Alessio et al. 1998). A point-like source equivalent to the star must present the same values for $J^*(0)$ and $H^*(0)$. The incidence angle is given by:

$$ \mu^* = H^*(0)/J^*(0). \quad (60) $$

One finally obtains, using the gray optical depth $\tau$, the same expression (33) of Malbet & Bertout (1991) except for the factor $\omega^*$ already discussed in the previous section,

$$ T_A^4 = \frac{(T_{\text{eff}})_A^4}{4f_K} \left[ \mu^* \omega^* \left(1 - e^{-\tau/\mu^*}\right) + \frac{f_K(0)}{f_H} + f_K \frac{e^{-\tau/\mu^*}}{\mu^*} \right] \quad (61) $$

where the effective temperature of the stellar reprocessing is:

$$ (T_{\text{eff}})_A^4 = 4\pi H^*(0)/\sigma_B. \quad (62) $$

### 3.3. Equation of state – Hydrostatic equilibrium

The equation of state in the disk is:

$$ P = f_{\text{gas}} c_s^2 \rho \quad (63) $$

$$ \text{with } c_s = \sqrt{\frac{k_B T}{M m_H}}. \quad (64) $$

First of all, radiation pressure is ignored in Eq. (63) as suggested by Shakura & Sunyaev (1973). If the disk is assumed to be homogeneous and mixed by the turbulence, the mass fraction of gas, $f_{\text{gas}}$, is a known function of $T$ and $\rho$, and might even be a constant if evaporation is ignored. The dust is only a few percent of the mass of the disk, so we use $f_{\text{gas}} = 1$ in the present work.

We also assume that the disk dynamics is dominated by the central star, so that it is Keplerian. Then, in the thin approximation, $G_z$ is a linear function of $z$. Combined with Eq. (63), the hydrostatic equilibrium and the local height scale become:

$$ \frac{d^2 P}{dm^2} = - \frac{c_s^2 \Omega^2}{P}, \quad (65) $$

$$ h = \frac{c_s}{\Omega}, \quad (66) $$

$$ \text{with } \Omega = \sqrt{\frac{GM_*}{r^3}}. \quad (67) $$

### 3.4. Self-similar accretion model

The radial distribution of energy dissipation for the accretion $U_V(r)$ and the surface density $\Sigma(r)$ are essential parameters for the initialisation step of our code. In order to be as general as possible, we use a prescription with a self-similar accretion rate, instead of the model by Shakura & Sunyaev (1973) which assumes a uniform accretion rate over the disk. We assume that some ejection process removes material away from the close environment of the star within a radius $r_0$. A more detailed description of this mechanism can be found in Ferreira & Pelletier (1995). The accretion rate is parametrized as:

$$ \dot{M}(r) = \begin{cases} \dot{M}_{\infty} \left(\frac{r}{r_0}\right)^{\xi} & \text{if } r < r_0 \\ \dot{M}_{\infty} & \text{if } r \ge r_0 \end{cases} \quad (68) $$

where $0 \le \xi \le 1$ is the ejection index (see Ferreira & Pelletier 1995) and $r_0$ the cut-off radius. When $\xi = 0$, this model corresponds exactly to the standard accretion model by Shakura & Sunyaev (1973).

The dissipation per area unit $U_V(r)$ and the mass column $\Sigma(r)$ can then be expressed as:

$$ U_V(r) = \frac{3GM_* \dot{M}(r)}{8\pi r^3} f(r) \quad (69) $$

$$ \Sigma(r) = \frac{\dot{M}(r)}{3\pi\nu_{\text{visc}}} f(r) \quad (70) $$

where $f(r)$ is a factor close to unity at large $r$. Using the same method as Shakura & Sunyaev (1973) we obtain:

$$ f(r) = \begin{cases} \frac{1}{1 + 2\xi} \left[ 1 - \left(\frac{r_*}{r}\right)^{\xi + \frac{1}{2}} \right] & \text{if } r < r_0 \\ 1 - \frac{2\xi}{1 + 2\xi} \sqrt{\frac{r_0}{r}} - \frac{1}{1 + 2\xi} \left(\frac{r_*}{r_0}\right)^{\xi} \sqrt{\frac{r_*}{r}} & \text{if } r \ge r_0. \end{cases} \quad (71) $$
<table><caption>**Table 1.** Parameters used in the model.</caption><tr><td>r<sub>*</sub> = 2 R<sub>☉</sub></td><td>M<sub>*</sub> = 0.5 M<sub>☉</sub></td></tr><tr><td>T<sub>*</sub> = 4000 K</td><td>T<sub>ISM</sub> = 15 K</td></tr><tr><td>r<sub>0</sub> = 100 r<sub>*</sub></td><td>α = 0.01</td></tr><tr><td>r<sub>min</sub> = 2.2 R<sub>☉</sub></td><td>A<sub>V</sub> = 1</td></tr><tr><td>d = 450 pc</td><td></td></tr></table>

### 3.5. The numerical code

Our set of equations are coupled and non-linear. They need to be solved numerically, within the framework presented below.

**The grid:** we use a logarithmic radial grid. For each radius we chose a non-uniform vertical grid so that the mass column coordinate $m(z)$ is logarithmically sampled.

**Initial conditions:** from the radial structure described in Sect. 2.1, we compute an isothermal vertical structure.

**Iterative method:** Knowing the temperature and mass column, the hydrostatic equilibrium is solved for. We deduce the density and pressure distributions. Then the temperature and optical depth are computed again using the development of Sect. 2.3.2. Subsequent quantities are derived, like the kinematic viscosity. The value of $M(r)$ is refined so that it matches Eq. (70) using the average kinematic viscosity. The mass column grid is computed again with the new value of $M(r)$. The latter step is performed as many times as necessary to reach convergence (typically in ten iterations).

The parameters used for the computations are reported in Table 1.

## 4. Results and discussion

In this paper, we chose to treat mainly the radiative transfer in the disk. Even if we have developed the formalism of the central source heating, we will ignore its contribution in a first step and concentrate only on the viscous dissipation. In a future paper, we will present the case of heating by the central star. Moreover the heating by the central star becomes dominant at radii larger than a few astronomical units as shown in Fig. 3 of D'Alessio et al. (1998) so the results presented here will be correct only for the inner part of the disk and therefore for the visible and infrared thermal domain.

First, we compute the physical conditions in the disk: temperature, density, optical depth. Then we use these physical conditions to compute astronomical observables like SEDs, millimetric images and interferometric complex visibilities and direct images.

### 4.1. Physical conditions in the disk

#### 4.1.1. Temperature and density

A map of density and temperature for a standard disk ($\xi = 0$) with an accretion rate of $10^{-5} M_{☉}/\text{yr}$ is displayed in Fig. 1. The shape of the temperature iso-contours illustrates that the disk temperature is almost constant in

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Iso-contour maps of a T Tauri disk with $\dot{M} = 10^{-5} M_{☉}/\text{yr}$. Dashed lines: lines of equal temperature; Dotted lines: lines of equal density; Solid line: limits of the optically thick region ($\tau = 1$); Gray zone: zone of instability with respect to the convection.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Radial distribution of the temperature in the central layer and of the effective temperature (located at $\tau = 1$).</figcaption></figure>

the optically thick part, whereas the temperature changes quickly with the altitude in the upper layer. As previously guessed, the disk appears roughly isothermal over the vertical dimension in its inner part.

Figure 2 displays the variation of the temperature with the radius. We verify that the effective temperature, i.e. at $\tau = 1$, is proportional to the dissipated flux: it follows
<figure><img src="image_5.png" /><figcaption>Fig. 3. Column density in an α-disk for two different accretion rates. Solid lines: $\xi = 0$; Dashed lines: $\xi = 0.5$.</figcaption></figure>

an $r^{-3/4}$ law. We stress that the central temperature is much higher than the effective temperature.

We check that the disk is optically thick even at large radii while the geometrical thickness of the disk remains much smaller than unity: if we define the disk surface as the region of optical depth unity, we find $z/r \approx 0.1-0.2$.

The geometry of the density contours displayed in Fig. 1 have a different shape depending whether they are located in the central optically thick part or in the upper optically thin layers. Figure 1 also displays the disk parts which are unstable to convection by computing the Schwarzschild criterion. For the time being, the consequences of this convection are not taken into account in our model.

### 4.1.2. Column density

Figure 3 shows the mass column density as a function of the radius in the case of a standard α-disk with two values of the ejection parameter $\xi$. We take the cut-off $r_0$ at 100 stellar radii, i.e. about 1 AU. Two values of the accretion rate have been considered: $10^{-6}$ and $10^{-4} M_{\odot}/\text{yr}$. As expected, the radial structure of $\Sigma(r)$ changes in the close environment of the star when $\xi \neq 0$.

Power-laws are often used by observers in the description of the radial structure of disks, though one can see in Fig. 4 that the column density of an α-disk cannot be described with such a law. However, on restricted domains of radius, the mass column in an α-disk can be described by the equation

$$ \Sigma(r) \approx \Sigma_i (r/1 \text{ AU})^{-q_i} \text{ kg m}^{-2} \text{ for } r_i \le r \le r_{i+1} \quad (72) $$

<figure><img src="image_6.png" /><figcaption>Fig. 4. Column density in a disk with different viscosity prescriptions. Solid lines: $\alpha = 10^{-3}, 10^{-2}$; Dotted lines: $\beta = 10^{-4}, 10^{-3}$.</figcaption></figure>

<table><caption>**Table 2.** Parameters computed from the model to approximate the surface density by a power law by pieces, sorted by decreasing order of radius domain.</caption><thead><tr><th colspan="3">M&#x0307; = 10<sup>-6</sup> M<sub>&#x2213;</sub>/yr</th><th colspan="3">M&#x0307; = 10<sup>-4</sup> M<sub>&#x2213;</sub>/yr</th></tr><tr><th>&Sigma;<sub>i</sub></th><th>q<sub>i</sub></th><th>r<sub>i</sub></th><th>&Sigma;<sub>i</sub></th><th>q<sub>i</sub></th><th>r<sub>i</sub></th></tr></thead><tbody><tr><td>1 &times; 10<sup>6</sup></td><td>1.5</td><td>40</td><td>&mdash;</td><td>&mdash;</td><td>&mdash;</td></tr><tr><td>4 &times; 10<sup>3</sup></td><td>0.0</td><td>9</td><td>2 &times; 10<sup>4</sup></td><td>0.0</td><td>67</td></tr><tr><td>1.5 &times; 10<sup>4</sup></td><td>0.6</td><td>1</td><td>2.5 &times; 10<sup>5</sup></td><td>0.6</td><td>5</td></tr><tr><td>1.5 &times; 10<sup>4</sup></td><td>1.2</td><td>0.11</td><td>9 &times; 10<sup>5</sup></td><td>1.4</td><td>0.75</td></tr><tr><td>4 &times; 10<sup>5</sup></td><td>-0.3</td><td>0.06</td><td>1.5 &times; 10<sup>6</sup></td><td>-0.5</td><td>0.45</td></tr></tbody></table>

where the coefficient $q_i$ ranges from -0.4 to 1.5 and $\Sigma_i$ is a constant. These coefficients are given in Table 2 for two values of the accretion rate.

We show in Fig. 4 the influence of $\alpha$ on the structure of the disk. These results are compared with the ones obtained with the $\beta$ viscosity prescription. The general behaviour is somewhat different except for the 1–10 AU region. Therefore observing inside 1 AU or outside 10 AU could be decisive in determining the viscosity law.

### 4.2. Astronomical observables

One of the aims of this work is to compute physical conditions around a young star so that we can calculate the electromagnetic field emerging from this region. This emission can then be analysed by several types of instruments: photometers, spectrographs, imaging cameras, interferometers, etc. The next sections give examples of astronomical observables that can be used to constrain our knowledge of the circumstellar environment.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Spectral energy distribution (λFλ vs. λ) for different disk models Solid lines: standard disk; Dashed lines: disk with ξ = 0.5; Dashed & dotted lines: standard disk with an inner hole of 8r*.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Spectral energy distribution (λFλ vs. λ) of the disk model as a function of inclination. Solid lines: pole-on disk. Dashed lines: inclination of i = 45°. Dotted lines: same inclination without the flaring effect (see text for details).</figcaption></figure>

### 4.2.1. Spectral energy distribution

The spectrophotometry technique gives access to the spectral energy distribution. We integrate along the line of sight the contribution of each disk layer to the emergent flux. In the examples taken for Fig. 5, the SED is significantly altered in the range 1–5 µm with the self-similar accretion model. The inner hot regions being depleted, the

<figure><img src="image_9.png" /><figcaption>Fig. 7. Emergent millimetric flux from a pole-on disk for several accretion rates and viscosity prescriptions. Upper curves: $\dot{M} = 10^{-4} M_{\odot}/yr$; Lower curves: $\dot{M} = 10^{-6} M_{\odot}/yr$; Solid lines: $\alpha = 10^{-2}$; Dashed lines: $\beta = 10^{-4}$.</figcaption></figure>

short wavelengths are less present. We assumed an overall extinction of $A_V = 1$ in the visible and an inner radius of the disk of $r_{min} = 1.1$ stellar radius. However, Fig. 5 shows that the SED is significatively different from the one expected from a standard model with a large inner gap.

Figure 6 shows the influence of disk flaring on the intermediate wavelength range of the SED. The difference between the two curves representing the SED of a 45°-inclined disk lies in the location of the photosphere, i.e. the location of the disk layer that emits most of the emergent flux. In the flaring case, the flux comes from the layer where $\tau = 1$, whereas in the so-called no flaring case, the flux comes from the disk equatorial plane but with the same temperature as the effective temperature. This latter case corresponds to the standard radial model where the vertical structure is not taken into account. Even without the reprocessing of the stellar photons, the effect of the vertical extension of the disk atmosphere is therefore not negligible at intermediate wavelengths.

### 4.2.2. Millimetric images

Radio observations can help us discriminate between the two viscosity prescriptions, and especially radio observations in the outer disk, since they are sensitive to the mass column. The disk being optically thin in the millimetric continuum, the flux is given by:

$$F_{\nu} \propto \int_{-\infty}^{+\infty} T(z) \rho(z) dz \approx \Sigma \bar{T}. \qquad (73)$$
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Visibilities from a disk model as a function of the baseline. Gray: $\dot{M} = 10^{-4} M_{\odot}/\text{yr}$; Black: $\dot{M} = 10^{-6} M_{\odot}/\text{yr}$. For a given baseline the visibility varies with the position angle; the domain between the minimum and maximum values is filled.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Closure phases as a function of the position angle for a triplet combination of the VLTI auxiliary telescopes. Solid lines: $\dot{M} = 10^{-4} M_{\odot}/\text{yr}$. Dotted line: $\dot{M} = 10^{-6} M_{\odot}/\text{yr}$ with the largest baseline.</figcaption></figure>

Figure 7 compares the millimetric emissivity for an α-disk and a β-disk with the same characteristics. At intermediate radii, the two emissivities are more or less equivalent, however at larger radii the emissivities can follow somewhat different power laws. $F_{\nu}$ is proportional to $r^{-3/2}$ for an α-disk over the full range, whereas a β-disk presents a variable power-law exponent. In the domain where the disk is thermalized with the ISM, i.e. at large radii, $F_{\nu} \propto r^{-1/2}$.

### 4.2.3. Interferometric observables

Probing the inner parts of disks requires high angular resolution techniques that can be achieved with present or soon-to-be operated interferometric instruments (Malbet & Bertout 1995). The more massive disks ($\dot{M} > 10^{-6} M_{\odot}/\text{yr}$) should be resolved by the VLTI and its 100m-baseline, but we should also detect non-zero closure phase.

Figure 8 shows the visibility of a disk as a function of the projected baseline. Since the image is not centro-symmetric, the visibility depends on the position angle. Figure 9 displays the closure phase as a function of the position angle for different telescope combinations in the VLTI for a disk with an important accretion rate ($\dot{M} = 10^{-4} M_{\odot}/\text{yr}$). As a comparison, the closure phase is displayed for a disk with moderate accretion rate with the widest combination of the VLTI auxiliary telescopes. Because of the natural flaring, a disk that is not seen pole-on does not appear centro-symmetric (as can be noticed in Fig. 10), therefore leading to a non-zero closure phase.

### 4.2.4. Direct images

Synthetic infrared images of the disk thermal emission are displayed in Fig. 10 with a field of 50 mas for a distance of $d = 450$ pc. As explained above, they are not centro-symmetric because of the flaring. So far, we have only access to interferometric observations at this scale, so they are only relevant as a step to produce interferometric observables.

## 5. Conclusion

We have presented a new method to model the vertical structure of accretion disks. This method leads to analytical formulae for the temperature distribution which help to understand the behaviour of the radiation propagated inside the disk. Our model includes two sources of energy dissipation: viscous heating and reprocessing of external radiation like that emitted by the central star. We have applied these analytical results to the case of T Tauri disks in a variety of conditions showing that the method is versatile:

- with a modified Shakura-Sunyaev prescription for a disk with a non-constant accretion rate (Ferreira & Pelletier 1995);

- with a viscosity prescription following either the standard $\alpha$ or the $\beta$ one proposed by Huré et al. (2001).

We are able to simulate the conditions of temperature and density in any part of the circumstellar environment and to compute astronomical observables like SED, optical and millimeter images or interferometric visibilities. This
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Synthetic images of the thermal emission at 2.2 µm for different inclinations of the disk (from 0° to 60° by steps of 10°). $\dot{M} = 10^{-6} M_{\odot}/yr$. The field is 50 mas.</figcaption></figure>

code could be used also as starting conditions in a Monte-Carlo multiple scattering code in order to derive polarization maps.

When completed, this code will offer to observers a tool based on physical parameters to interpret their observations better than ad-hoc models based on power laws often used today.

*Acknowledgements.* We would like to thank E. di Folco who has partly worked on the code. We are grateful to C. Bertout, J. Bouvier, C. Ceccarelli, C. Dougados and F. Ménard for helpful discussions. We also thank the referee, Dr. Hubeny, for useful suggestions.

## Appendix A: Vertical structure of the temperature with several energy sources

Each source of energy locally contributes to the heating $u_k$ (see Sect. 2.3.4). The total local heating is then $u = \sum_k u_k$ and the energy dissipated per unit of disk surface is then $U = \sum_k U_k$ with

$$U_k M = \int_0^M u_k(\zeta) \, d\zeta. \qquad (A.1)$$

We now write the vertical distribution of energy dissipation, $\theta$ as a function of the different vertical distributions of the energy source, $\theta_k$:

$$\begin{aligned} \theta(m) &= \int_0^m u(\zeta) \, d\zeta / (UM) \\ &= \sum_k \frac{U_k}{U} \theta_k(m) \end{aligned} \qquad (A.2)$$

with $\theta_k(m) = \int_0^m u_k(\zeta) \, d\zeta / (U_k M)$.

Similarly for $\Delta\tau_H$:

$$\begin{aligned} \Delta\tau_H(m) &= \int_0^m \theta(\zeta)\tau_H(\zeta) \, d\zeta \\ &= \sum_k \frac{U_k}{U} (\Delta\tau_H)_k \end{aligned} \qquad (A.3)$$

with $(\Delta\tau_H)_k(m) = \int_0^m \theta_k(\zeta)\tau_H(\zeta) \, d\zeta$.

Since,

$$\frac{u(m)}{U} = \sum_k \frac{U_k}{U} \frac{u_k(m)}{U_k} \qquad (A.4)$$

and $\sum_k (U_k/U) = 1$, we can rewrite Eq. (26) by replacing each term in the brackets by a sum over the $k$ sources:

$$T^4 = \frac{\kappa_J T_{\text{eff}}^4}{4\kappa_B f_K} \sum_k \frac{U_k}{U} \left[ \left( \tau_H - (\Delta\tau_H)_k + \frac{f_K(0)}{f_H} \right) + \frac{f_K}{M\kappa_J} \frac{u_k}{U_k} \right]. \qquad (A.5)$$

If we define the effective temperature of the $k$th energy source,

$$\sigma_B (T_{\text{eff}})_k^4 = U_k M, \qquad (A.6)$$

then Eq. (26) can be written as the following sum:

$$T^4 = \sum_k \frac{\kappa_J (T_{\text{eff}})_k^4}{4\kappa_B f_K} \left[ \left( \tau_H - (\Delta\tau_H)_k + \frac{f_K(0)}{f_H} \right) + \frac{f_K}{M\kappa_J} \frac{u_k}{U_k} \right]. \qquad (A.7)$$

Therefore the vertical distribution of the temperature for a disk with multiple energy dissipation sources is the same as the sum of the individual temperature distributions using the global radiation in the disk, i.e. with $J$, $H$, $K$, $\tau_H$ computed from the global temperature distribution.

## Appendix B: Incoming radiation field

We show in this appendix how to compute $J^0(m)$, $J^s(m)$, $H^0(m)$, and $H^s(m)$ in the case of an external point-like source. The intensity of the source is $J^0$, giving us a boundary condition to solve the equations of transfer (44) and (45). If we define the optical depths,

$$\tau^0 = \int_0^m \chi_J^0(\zeta) \, d\zeta \qquad (B.1)$$

$$\tau^s = \int_0^m \kappa_J^s(\zeta) \, d\zeta \qquad (B.2)$$

we obtain:

$$J^0(m) = J(0)e^{-\tau^0(m)/\mu^0} \qquad (B.3)$$

$$J^s(m) = \Lambda_{\tau^s(m)} \left( \frac{\sigma^0}{\kappa^s} J^0 \right) \qquad (B.4)$$

where $\Lambda$ is a commonly used operator (cf. Mihalas 1978) defined by:

$$\Lambda_{\tau}(f) = \frac{1}{2} \int_{0}^{+\infty} E_{1} |t - \tau| f(t) \, dt \qquad (B.5)$$

$$\text{and } E_n(t) = \int_{1}^{+\infty} u^{-n} e^{-ut} \, du. \qquad (B.6)$$
The terms $H^0(m)$ and $H^s(m)$ can be computed in the same way as $J^0(m)$ and $J^s(m)$:

$$H^0(m) = \mu^0 J(0)e^{-\tau^0(m)/\mu^0} \qquad (B.7)$$

$$H^s(m) = \Phi_{\tau^s(m)} \left( \frac{\sigma_J^0}{\kappa_J^s} J^0 \right) \qquad (B.8)$$

with

$$\Phi_{\tau}(f) = 2 \int_{0}^{+\infty} E_2|t - \tau|f(t) \operatorname{sgn}(t - \tau) dt \qquad (B.9)$$

where sgn($x$) is the sign function.

## References

* Adams, F. C., Shu, F. H., & Lada, C. J. 1988, ApJ, 326, 865

* Basri, G., & Bertout, C. 1989, ApJ, 341, 340

* Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987

* Bertout, C., Basri, G., & Bouvier, J. 1988, ApJ, 330, 350

* Burrows, C. J., Stapelfeldt, K. R., Watson, A. M., et al. 1996, ApJ, 473, 437

* Calvet, N., Patio, A., Magris, G., & D'Alessio, P. 1991, ApJ, 380, 617

* Casse, F., & Ferreira, J. 2000a, A&A, 353, 1115

* Casse, F., & Ferreira, J. 2000b, A&A, 361, 1178

* Chandrasekhar, S. 1960, in Radiative Transfer (New York: Dover)

* Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368

* Chiang, E. I., & Goldreich, P. 1999, ApJ, 519, 279

* D'Alessio, P., Canto, J., Calvet, N., et al. 1998, ApJ, 500, 411

* D'Alessio, P., Calvet, N., Hartmann, L., et al. 1999, ApJ, 527, 893

* Di Folco, E. 1999, graduation dissertation from Grenoble University

* Feautrier, P. 1964, C. R. Acad. Sci. Paris, 258, 3189

* Ferreira, J., & Pelletier, G. 1995, A&A, 295, 807

* Franck, J., King, A. R., & Raine, D. J. 1985, in Accretion Power in Astrophysics (Cambridge: University Press)

* Hubeny, I. 1990, ApJ, 350, 632

* Hubeny, I. 1991, The 6th Institute d'Astrophysique de Paris (IAP) Meeting: Structure and Emission Properties of Accretion Disks, IAU Colloq., 129, 227

* Huré, J.-M., Richard, D., & Zahn, J.-P. 2001, A&A, 367, 1087

* Lazareff, B., Monin, J.-L., & Pudritz, R. E. 1990, ApJ, 358, 170

* Malbet, F., Lachaume, R., Monin, J.-L., & Berger, J.-P. 2000, SPIE, 4006, 243

* Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

* Malbet, F. 1992, Ph.D. from Paris University

* Malbet, F., & Bertout, C. 1995, A&AS, 113, 369

* Malbet, F., & Bertout, C. 1991, ApJ, 383, 814

* Meyer, F., & Meyer-Hofmeister, E. 1982, A&A, 106, 34

* Mihalas, D. 1978, in Stellar atmospheres (San Francisco: Freeman)

* Monin, J.-M., & Bouvier, J. 2000, A&A, 356, L75

* Pringle, J. E. 1981, ARA&A, 19, 137

* Ruden, S. P., & Pollack, J. B. 1991, ApJ, 375, 740

* Rydgren, A. E., & Zak, D. S. 1987, PASP, 99, 191

* Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

* Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781

* Stapelfeldt, K. R., Krist, J. E., Ménard, F., et al. 1998, ApJ, 502, L65

* Terquem, C., Papaloizou, J. C. B., Nelson, R. P., & Lin, D. N. C. 1998, ApJ, 502, 788

* Terquem, C. 1998, ApJ, 509, 819

* Whitney, B. A., & Hartmann, L. 1992, ApJ, 395, 529