3 Radiative transfer

In this Section the continuum and molecular line radiative transfer codes used to constrain the physical and chemical structure of the envelope around IRAS 16293-2422 are presented. The envelope is assumed to be spherically symmetric and the approach adopted here is to first determine the density and temperature structures from the dust modelling. This, in turn, allows for abundances of various molecules present in the circumstellar envelope to be determined.

   
3.1 Dust radiative transfer model

In order to model the observed continuum emission and to be able to extract some basic parameters of the dusty envelope around IRAS 16293-2422 the publically available dust radiative transfer code DUSTY (Ivezic et al. 1999) has been adopted. DUSTY makes use of the fact that in some very general circumstances the radiative transfer problem of the dust possesses scaling properties (Ivezic & Elitzur 1997). The solution is presented in terms of the distance, r/ $r_{{\rm i}}$, scaled with respect to the inner boundary $r_{{\rm i}}$. In addition to the properties of the dust and the relative size of the envelope, the only parameter needed for a full description of the problem is the spectral shape of the radiation emitted by the central source. This means that the luminosity is totally decoupled from the radiative transfer problem and it is only used to scale the solution in order to obtain the absolute distance scale. This scaling property of the dust radiative transfer is very practical, in particular when modelling a large number of similar sources, and is put to use in Jørgensen et al. (2002) for a survey of low mass protostars at various evolutionary stages and in Hatchell et al. (2000) to model a selection of high mass protostars.

The most important parameter controlling the output is the dust optical depth

$$\tau_{\lambda} = \kappa_{\lambda} \int_{r_{\rm i}}^{r_{\rm e}} \rho_{\rm d}(r)~{\rm d}r,$$ (1)

where $\kappa_{\lambda}$ is the dust opacity, $\rho_{\rm d}$ the density distribution of the dust in the envelope, and $r_{{\rm i}}$and $r_{\rm e}$ the inner and outer radius of the dust envelope, respectively. Introducing the dust-to-gas mass ratio, $\delta$, allows Eq. (1) to be written

$$\tau_{\lambda} = \kappa_{\lambda} \delta m_{\rm H_2} \int_{... ...\rm d}r = \kappa_{\lambda} \delta {\rm m_{H_2}} N_{\rm H_2},$$ (2)

where $n_{\rm H_2}$ is the number density distribution of molecular hydrogen, $m_{\rm H_2}$ the mass of an hydrogen molecule, and $N_{\rm H_2}$ the column density of molecular hydrogen. In the derivation of Eq. (2) any possible drift velocity between the dust and the gas has been neglected. In what follows $\delta=0.01$ is assumed. The dust opacities from Ossenkopf & Henning (1994) were used, corresponding to coagulated dust grains with thin ice mantles, at a density of $n_{\rm H_2} \sim 10^6$ cm^-3 (Col. 5 in their Table 1, hereafter OH5). Van der Tak et al. (1999) considered various sets of dust optical properties when modelling the high mass young stellar object GL 2591, and found that models using the OH5 opacities were the only ones that gave envelope masses consistent with those derived from the modelling of the molecular line emission. In Sect. 4.2 another set of opacities will also be considered, more appropriate for regions where the ices have evaporated from the dust grains.

The dust temperature at the inner radius of the envelope is fixed to 300 K and this sets the inner radius. The choice of this temperature is motivated by the observations of line emission arising from highly excited molecules in the envelope. However, the inner regions are complex with breakdown of spherical symmetry and interactions between disk, envelope and outflow. In the present analysis these complications are ignored and for simplicity a smooth and spherically-symmetric envelope is assumed. Recently, Ceccarelli et al. (2000a) successfully modelled molecular line emission in the envelope around IRAS 16293-2422 down to $\sim$30 AU, assuming spherical symmetry. While the separation of the two protostars is approximately 800 AU (Looney et al. 2000), one of the protostars IRAS 16293A (MM1) exhibits jet-like centimetre wavelength emission, water maser emission and associated millimetre molecular emission, thus appearing to be significantly more active (Wootten & Loren 1987; Mundy et al. 1992; Schöier et al. 2002a, in prep.), which justifies the approach taken here. The central source of radiation is assumed to arise from a blackbody at 5000 K. This is an oversimplification considering the binary nature of IRAS 16293-2422 and the uncertainty in the intrinsic SED of a protostar. The final model does not depend on the exact stellar temperature adopted, however, within a reasonable range of values, since the radiation is totally reprocessed by the circumstellar dust. The input parameters are summarized in Table 4.

The observational constraints, as presented in Sect. 2, are the SED and radial brightness distributions at 450 $\mu$m and, to a lesser extent, 850 $\mu$m. The ability of the model to reproduce the observational constraints are quantified using the chi-squared statistic

$$\chi^2 = \sum^N_{i=1} \left [ \frac{(F_{{\rm mod},i}-F_{{\rm obs},i})}{\sigma_i}\right ]^2,$$ (3)

where F is the flux and $\sigma_i$ the uncertainty in observation i, and the summation is done over all N independent observations. The radial brightness distributions from DUSTY are extended into 2D surface brightness maps and convolved with the beam as determined from planet observations. The beam convolved maps are then azimuthally averaged in 3$\arcsec$ bins. In the $\chi ^2$ fitting procedure only data points separated by one full beam are used since the $\chi ^2$-analysis, in practice, requires uncorrelated measurements. In the analysis only the inner 50$\arcsec$ of the brightness distributions are considered which are well above the background emission and should not be significantly affected by the 120$\arcsec$ chop throw. Furthermore, the IRAS fluxes from the model were convolved with the proper filters before the $\chi ^2$ analysis of the SED was made. The results from the dust radiative transfer are presented in Sect. 4.

  

Table 4: Summary of the dust radiative transfer analysis of IRAS 16293-2422 using a single power-law density distribution (see text for details).

   
3.2 Molecular line radiative transfer model

In order to derive accurate molecular abundances for the wealth of molecular line emission detected toward this source the detailed non-LTE radiative transfer code of Schöier (2000), based on the Monte Carlo method, was used. The code produces output that is in excellent agreement with the Monte Carlo code presented by Hogerheijde & van der Tak (2000). It has also been tested against other molecular line radiative transfer codes, for a number of benchmark problems, to a high accuracy (van Zadelhoff et al. 2002).

Adopting the parameters of the circumstellar envelope derived from the dust radiative transfer analysis, the Monte Carlo code calculates the steady-state level populations of the molecule under study, using the statistical equilibrium equations. In the Monte Carlo method, information on the radiation field is obtained by simulating the line photons using a number of model photons, each representing a large number of real photons from all transitions considered. These model photons, emitted locally in the gas as well as injected from the boundaries of the envelope, are followed through the envelope and the number of absorptions are calculated and stored. Photons are spontaneously emitted in the gas with complete angular and frequency redistribution, i.e., the local emission is assumed to be isotropic and the scatterings are assumed to be incoherent. The weight of a model photon is continuously modified as it travels through the envelope, to take the absorptions and stimulated emissions into account. When all model photons are absorbed in, or have escaped from, the envelope the statistical equilibrium equations are solved and the whole process is then repeated until some criterion for convergence is fulfilled. Once the molecular excitation, i.e., the level populations, is obtained the radiative transfer equation can be solved exactly. The resulting brightness distribution is then convolved with the appropriate beam to allow a direct comparison with observations. In this analysis, the kinetic temperature of the gas is assumed to follow that of the dust (Ceccarelli et al. 1996; Doty & Neufeld 1997; Ceccarelli et al. 2000a).

Typically, energy levels up to $\sim$500 K in the ground vibrational state were retained in the analysis. Vibrationally excited levels are not included since the radiative excitation due to the dust is generally inefficient and, for the temperature and density ranges present here, collisional excitation to these levels is negligible. Line emission from rotational transitions within vibrationally excited states has, however, been observed for some species (e.g., the CS (v =1, $J =7\rightarrow 6$) Blake et al. 1994). In Sect. 5 the nature of such emission is discussed further. Collisional rate coefficients are taken from the literature and, in the case of some linear molecules, extrapolated both in temperature and to transitions involving energy levels with higher J quantum numbers when needed (see Schöier et al. 2002b, in prep. for details). For other (non-linear) species for which such extrapolations are not obvious, the excitation is assumed to be in LTE for levels for which no collisional rate coefficients are available.

In what follows, all quoted abundances, $f_{\rm X}$, for a particular molecular species X are relative to that of molecular hydrogen, i.e.,

$$f_{\rm X}(r) = \frac{n_{\rm X}(r)}{n_{{\rm H_2}}(r)},$$ (4)

and are initially assumed to be constant throughout the envelope. Subsequently, it will be shown that in order to model the observed line emission for some molecules, e.g., H₂CO and SiO, this latter constraint has to be relaxed and a jump in f needs to be introduced (see also Ceccarelli et al. 2000a,b). Van der Tak et al. (2000b) used a similar approach for the case of massive protostars.

The beam profile used in the convolution of the modelled emission is assumed to be Gaussian which is appropriate at the frequencies used here. The best fit model is estimated from the $\chi ^2$-statistic defined in Eq. (3) using the observed integrated intensities and assuming a 30% calibration uncertainty. Even though the new data presented in Table 3 appear to be slightly better calibrated, the old data set constitutes the vast majority of observational constraints so the larger calibration uncertainty is applied to the full data set, for simplicity. Only in the cases of CO and CS was a lower calibration uncertainty of 15% adopted. These molecules are regularly observed towards IRAS 16293-2422 and used as standard spectra.