A&A 418, 1021-1034 (2004)
DOI: 10.1051/0004-6361:20034476

Physical-chemical modeling of the low-mass protostar IRAS 16293-2422

S. D. Doty ¹ - F. L. Schöier ^2,3 - E. F. van Dishoeck ³

1 - Department of Physics and Astronomy, Denison University, Granville, OH 43023, USA
2 - Stockholm Observatory, AlbaNova, 10691 Stockholm, Sweden
3 - Leiden Observatory, PO Box 9513, 2300 RA Leiden, The Netherlands

Received 8 October 2003 / Accepted 29 January 2004

Abstract
We present detailed gas-phase chemical models for the envelope of the low-mass star-forming region IRAS 16293-2422. By considering both time- and space-dependent chemistry, these models are used to study both the physical structure proposed by Schöier et al. (2002), as well as the chemical evolution of this region. A new feature of our study is the use of a detailed, self-consistent radiative transfer model to translate the model abundances into line strengths and compare them directly with observations of a total of 76 transitions for 18 chemical species, and their isotopes. The model can reproduce many of the line strengths observed within 50%. The best fit is for times in the range of $3\times 10^3{-}3\times 10^4$ yrs and requires only minor modifications to our model for the high-mass star-forming region AFGL 2591. The ionization rate for the source may be higher than previously expected - either due to an enhanced cosmic-ray ionization rate, or, more probably, to the presence of X-ray induced ionization from the center. A significant fraction of the CO is found to desorb in the temperature range of 15-40 K; below this temperature $\sim$90% or more of the CO is frozen out. The inability of the model to explain the HCS⁺, C₂H, and OCS abundances suggests the importance of further laboratory studies of basic reaction rates. Finally, predictions of the abundances and spatial distributions of other species which could be observed by future facilities (e.g. Herschel-HIFI, SOFIA, millimeter arrays) are provided.

Key words: stars: formation - stars: individual: IRAS 16293-2422 - ISM: molecules

1 Introduction

The distribution and composition of dust and gas around isolated low-mass young stellar objects (YSOs) is receiving increased attention both observationally and theoretically. While the general process of low mass star formation is relatively well understood (see e.g. Shu et al. 1993; Evans 1999, and others), many details on the chemical and physical structure at different stages of evolution remain uncertain. In particular, the warm and dense gas in the very interior of these star-forming regions provides a rich chemical environment with which to probe their structure, properties, and evolution. Beyond their own intrinsic interest, these regions may provide a link to the so-called hot cores observed toward many high-mass star-forming regions (e.g., Walmsley & Schilke 1993).

Rapid advances in both observational and modeling capabilities allow much more quantitative studies of the chemistry in YSO envelopes than was possible even a few years ago. Several different steps can be distinguished (see Fig. 1). Thanks to the advent of large-format bolometer arrays, most studies nowadays start with an analysis of the spatial distribution of the submillimeter continuum emission from dust and its spectral energy distribution (SED) (e.g., Shirley et al. 2000; Jørgensen et al. 2002; Schöier et al. 2002). Through continuum radiative transfer calculations, both the density profile  $n\propto r^{-p}$ as a function of radius r and the dust temperature structure  $T_{\rm dust}(r)$ can be determined self-consistently. For the gas, two approaches can subsequently be taken. In atmospheric chemistry, these two cases are commonly known as the "forward'' and "backward'' or "retrieval'' methods. In the "empirical model'' ("retrieval''), $T_{\rm gas}$ is taken to be equal to  $T_{\rm dust}$and the excitation, radiative transfer and fluxes of the various molecular lines are calculated for an assumed abundance profile  $x(r)=n({\rm X})(r)/n$(H₂)(r). This trial abundance profile is then varied until the best agreement with observations is obtained. In practice, only two types of abundance profiles are considered: a constant abundance throughout the envelope or a "jump'' profile in which the abundance is increased by a large factor in the inner warm region due to ice evaporation. Such models have successfully been applied to both high- (e.g., van der Tak et al. 2000) and low-mass YSOs (e.g., Ceccarelli et al. 2000b; Schöier et al. 2002), and work best for molecules for which a large set of lines originating from levels with a range in energy has been observed. This method provides abundances for comparison with chemical models, but does not test the chemical networks directly.

Figure 1: An overview of methods used for constraining the physical and chemical structure of YSO envelopes. In general, the steps include the determination of the physical structure through dust modeling, calculation of gas temperature, and adoption of a chemical model. This combination produces observables (column densities, line fluxes, etc.) to compare with observations. The source parameters are determined by adjusting them until a best fit is obtained. Two important points of divergence include the use of self-consistent vs. approximate radiative transfer, and the use of a time-dependent chemical network ("Full chemical model'') vs. simple trial abundances ("Empirical model'') (figure adapted from van Dishoeck & van der Tak 2000, and van Dishoeck 2003).

Open with DEXTER

The second, ab initio or "forward'' approach is the "full chemical model'', in which only the density structure derived from the dust is adopted as a starting point. Given initial abundances and a detailed chemical network, the abundances of various molecules can be solved as functions of position and time and the gas temperature can be calculated explicitly by solving the full thermal balance of the gas. The physical structure n(r) and T(r) can either be taken to be constant with time or to vary according to some (dynamical) prescription. Such time- and space-dependent chemical models have been applied to low-mass YSOs by Ceccarelli et al. (1996, hereafter CHT96) and Rodgers & Charnley (2003), and to specific high-mass sources by Millar et al. (1997; G34.3+0.15) and Doty et al. (2002; AFGL 2591). The main quantities to be determined are the best-fit time (or, in a dynamical model, mass-accretion rate) and other parameters which enter the chemical models, such as the cosmic ray ionization rate. The chemical models themselves can be tested by comparing the abundance profiles  $x(r,t_{\rm fit})$ at the best-fit time with those derived through the empirical method. In this way, they also provide a guideline for more complicated abundance profiles to adopt in the empirical method.

In this paper, we describe a detailed "full chemical model'' of the best studied low-mass YSO, IRAS 16293-2422. A novel feature is the addition of a full Monte Carlo radiative transfer calculation of the resulting line fluxes for direct comparison with observations (see bottom-right part of Fig. 1). Such models provide the most complete test of our understanding of the physical and chemical structure of YSO envelopes.

IRAS 16293-2422 is a nearby ($\sim$160 pc, Whittet 1974) low mass, low luminosity (27 $L_{\odot}$) protostellar object located within the $\rho$ Ophiuchus molcular cloud. It has an exceptionally rich and well-studied spectrum (e.g., Blake et al. 1994; van Dishoeck et al. 1995; Cecarelli et al. 2000a,b; Schöier et al. 2002), and is therefore considered the prototypical low-mass source for chemical studies, much like Orion is for high-mass objects. Ceccarelli et al. (2000a,b) used the physical structure based on CHT96 combined with a restricted chemical network to analyze data of H₂CO, H₂O, and SiO in a "full chemical model'', and found strong evidence for large abundance enhancements of these species in the innermost part ($\leq$150 AU) of the envelope. The evaporated species may subsequently drive a complex "hot core'' chemistry leading to the even more complex organic molecules which have recently been detected in IRAS 16293-2422 (Cazaux et al. 2003).

In a later analysis, Schöier et al. (2002) combined dust/SED modeling of the physical structure, multiple line observations covering a range of excitation conditions, and a detailed radiative transfer analysis in an "empirical model'' to infer the structure of IRAS 16293-2422. This work supported the general conclusion of a "hot core'' where the abundances of key molecules are enhanced by several orders of magnitude due to evaporation of ices. The model employed only uniform and "jump'' abundances, however, which may not be representative of the detailed time- and space-dependent chemistry.

Recently, Doty et al. (2002) described such a time- and space-dependent physical/chemical model for static YSO envelopes including the hot core chemistry. By combining the model results with observations of many species of one particular high-mass YSO, AFGL 2591, it has been shown that it may be possible to not only confirm the gross source structure, but also constrain source properties such as age, ionization rate, and role of grains in determining the chemical structure (see also Boonman et al. 2003). Here the "full chemical model'' of Fig. 1 was adopted, but the self-consistent line radiative transfer was performed for only a subset of the species.

In this paper, we report on the application of the physical/chemical model of Doty et al. (2002) to the low-mass YSO IRAS 16293-2422. These results are combined with a self-consistent radiative transfer model, and applied to the full multi-species, multi-transition dataset of Schöier et al. (2002). By comparison with the case of AFGL 2591, we can also directly determine the differences in derived model parameters for a low- and a high-mass YSO (van Dishoeck 2003). The models and observations are briefly described in Sect. 2. The observations are then used with the models to constrain the source properties in Sect. 3. Finally, we summarize the results and conclude in Sect. 4.

2 Existing observations and models

2.1 Observations

IRAS 16293-2422 has been well-observed both in the continuum and in various submillimeter molecular lines. While no new observations are presented in this paper, we briefly note and discuss the observational data as they provide the constraints placed on the model.

The SED of IRAS 16293-2422 in the range $60{-}2900~{\mu}$m is presented by Schöier et al. (2002). High angular resolution JCMT data allowed the determination of radial brightness distributions at 450 and 850 ${\mu}$m on scales of 9'' and 15'' (1400 and 2400 AU) respectively.

Table 1: IRAS 16293-2422 physical structure from dust modeling by Schöier et al. (2002).

The molecular line data utilized here are taken primarily from the large surveys of IRAS 16293-2422 by Blake et al. (1994) and van Dishoeck et al. (1995). Additional complementary data are taken from the JCMT public archive (Schöier et al. 2002). These data are supplemented by the H₂CO lines of Loinard et al. (2000). The data set - a total of 76 transitions for 18 species considered here - has the advantage that it samples the full radial range of the envelope, providing probes over a wide range of physical, thermal, and chemical conditions. Only information on the lowest transitions of the molecules, which occur at millimeter wavelengths and probe the very coldest outer parts and surrounding cloud, is lacking. In general, the calibration uncertainty of each individual line is $\sim$30%.

2.2 Model

Here a brief synopsis of the physical, thermal, chemical, and radiative transfer models are provided. For more detailed information, see Doty et al. (2002), Schöier et al. (2002), Doty & Neufeld (1997), and references there.

2.2.1 Physical and thermal structure

We adopt a spherically symmetric static model of the extended envelope of IRAS 16293-2422 surrounding the two protostellar sources in the center (Looney et al. 2000). The observational continuum data were combined and simultaneously modeled by Schöier et al. (2002) with the publicly available radiative transfer code DUSTY (Ivezic et al. 1999). This analysis allows the source structure properties (e.g., envelope size, density power law, continuum optical depth) to be determined to within approximately $\pm$20% (see also Doty & Palotti 2002). The adopted source properties are presented in Table 1, and the density and temperature structure are reproduced in Fig. 2.

The successful line modeling of Ceccarelli et al. (2000a) and Schöier et al. (2002) in spherical symmetry down to $\sim$30 AU suggests that the assumption of spherical symmetry may be largely justifiable. Most chemical data were obtained with a 15'' beam, which yields a linear size of 2400 AU at the assumed distance of 160 pc, larger than the 800 AU separation of the central protostars. Schöier et al. (2003) recently used interferometric observations to determine the structure below 1000 AU, finding that the binary has cleared out most of the material in the inner part of the envelope, and that there exist two unresolved central sources with best-fit disk sizes of $\sim$250 AU in diameter. However, they find that the while the detailed results for the inner envelope leave the inner (r < 400 AU) structure somewhat uncertain, their results have little effect on the extended envelope (r > 400 AU).

	Figure 2: Physical and thermal structure of IRAS 16293-2422. The density and temperature are from the model of Schöier et al. (2002). The gas temperature is assumed to follow the dust temperature.
Open with DEXTER

A majority of the modeled envelope exists at relatively high densities and optical depths, leading to a strong thermal coupling between the gas and dust. As a result, the gas temperature is assumed to follow the dust temperature (Ceccarelli et al. 1996; Doty & Neufeld 1997; Ceccarelli et al. 2000a). Test calculations have shown that this assumption is sufficient for both the chemistry (Doty et al. 2002) and radiative transfer through molecular lines (Boonman et al. 2003).

2.2.2 Chemistry

The chemical model is based upon the UMIST gas-phase chemical reaction network (Millar et al. 1997, hereafter MFW), including reactions to model the hot-core chemistry. Pseudo time-dependent models of the chemical evolution over 30 radial grid points were constructed, providing a time- and space-dependent chemical evolution. The local parameters (hydrogen density, temperature, and optical depth) at each radial point are taken from the physical and thermal structure calculations above. For the majority of species, the initial abundances of the high-mass source AFGL 2591 (Doty et al. 2002), are utilized as shown in Table 2. The chemistry of deuterated molecules is not considered.

The effects of freeze out onto and desorption from dust grains are approximated. Instead of explicit freeze out or desorption with time, the desorption is taken to be instantaneous at 100 K, where expected grain mantle species are injected into the gas (Charnley 1997), in keeping with the timescales observed in the laboratory (Fraser et al. 2001). Below 100 K, we deplete the gas phase abundance of many species expected to be in ices (e.g., H₂O). The only exceptions are CO, which we take to desorb at  $T_{{\rm CO}}$, and H₂CO and CH₃OH, which we take to desorb or "jump'' at T_X.

The effects of photodissociation from the interstellar radiation field at the outer boundary is included, but is generally small due to the high optical depth, and the coarseness of the spatial grid considered. Finally, we have considered the effect of metal depletion by significantly reducing the initial Fe abundance. We find that metal depletion makes only a small difference, worsening the fits by only a few percent on average.

Table 2: Initial abundances at t=0 relative to H₂ for AFGL 2591.

2.2.3 Radiative transfer

The molecular line radiative transfer is accomplished through a non-LTE, Monte-Carlo model described in Schöier (2000). This code has been benchmarked to high accuracy against a wide range of other molecular line radiative transfer models (van Zadelhoff et al. 2002). In this model, the spatial molecular abundances x(r,t) are combined with the adopted physical structure to compute the excitation and resulting line profiles for all transitions up to $\sim$500 K in the ground vibrational state of the observed molecules. Chemical evolution times from $3 \times 10^{2}$ years to $3 \times 10^{5}$ years are considered, with one dex spacing.

3 Results

In this section the results of our physical/thermal/ chemical modeling of IRAS 16293-2422 are presented, and the comparison of line strengths predicted from this model to those observed. As a metric of the comparison, we adopt the mean percentage magnitude difference between the predicted and observed line strengths, given by

 

$$\Delta = \frac{1}{N_{\rm {lines}}} \sum_{i=1}^{N_{{\rm lines}}} \left\vert \frac{F_{{\rm mod},i}-F_{{\rm obs},i}}{F_{{\rm obs},i}} \right\vert\cdot$$ (1)

This form has the advantage that it measures the size of the difference, without allowing equally balanced high and low values to cancel out. Furthermore, the summation over lines implicitly provides a greater weight to molecules for which more data - and thus more constraints - exist. Measures using the sum over molecules instead of lines show qualitatively similar, though often accentuated, results. We note that we adopt this measure instead of a $\chi^2$ analysis, as agreement to within a factor of 10 are considered good for chemical modeling. This value is enough larger than the statistical uncertainty in the line observations so as to invalidate the statistical meaning of $\chi^2$ measure.

The parameters varied in the models are the cosmic ray ionization rate $\zeta$, the adopted initial abundances in the inner and outer regions, and the desorption temperatures of selected species (CO, H₂CO, CH₃OH). Detailed radial profiles of selected species are presented in Sects. 3.7 and 3.8.

3.1 General results

As discussed above, our base model is taken after the high-mass hot core + envelope of AFGL 2591 by Doty et al. (2002). A parameter space search, guided by results from previous studies, suggests a best fit to the observed data of IRAS 16293-2422 with only minor modifications to the AFGL 2591 initial chemical conditions.

Figure 3: The ratio of the predicted and observed line strengths for molecules observed toward IRAS 16293-2422. The errorbars represent the range of values between $3 \times 10^{3} < t({\rm yrs}) < 3 \times 10^{4}$. This represents the best fit time range, and is consistent with previous values based upon estimates of the infall rate and central mass. Guidelines at factor of 3 (dotted) and 10 (dashed) ratios are given to indicate good and acceptable fits.

Open with DEXTER

A comparison between the best fit model and observations is shown in Fig. 3. Here the ratio of the predicted to observed line strengths for each molecule observed is plotted. We find a best-fit time of $3 \times 10^{3} < t({\rm yrs}) < 3 \times 10^{4}$, with these times forming the range in the figure. These times are consistent with the age inferred by Schöier et al. (2002) from fitting a collapse model to the line profiles, and by using the constant infall rate of Ceccarelli et al. (2000a) with their preferred central mass of  $0.8~M_{\odot}$.

As can be seen in the figure, the majority of species (11 of 18) are fit to within 50% of the observations, thirteen are fit to within a factor of three, and 15 are fit to within a factor of 10, a level usually considered acceptable agreement in chemical modeling (see e.g., Millar & Freeman 1984; Brown & Charnley 1990; Terzieva & Herbst 1998). An interesting case is the $\sim$$20{-}30\%$ deviation in ¹³CO, and the small uncertainties on ¹³CO, C¹⁸O, and HCO⁺. It is possible that the ¹³CO discrepancy could be due to a different ¹²C/¹³C ratio than taken here. It may also be due to deviation of the real structure from the continuous, spherical symmetry we have adopted. In any case, the deviation is no larger than the expected calibration uncertainty of $\sim$$20{-}30\%$.

In the comparison there exist three species which deviate by more than an order of magnitude. The outliers are: OCS, C₂H, and HCS⁺, which have individual deviations of a factor 30, 30, and 100 respectively. While there is variation, our model tends toward producing too little emission. The difference between the model and observations is  $\Delta=0.51$ when these are omitted. Including them raises $\Delta$ to 0.87. This quantitatively confirms the agreement between the model and observations at the level of a factor of two for most species.

Table 3: Differences between best fit model for IRAS 16293-2422 and AFGL 2591.

The differences between the AFGL 2591 model and the best fit model here are summarized in Table 3. As can be seen, the differences are generally minor. They are: (1) a large increase in the cosmic-ray ionization rate by a factor of >10 over the "standard'' value of $\sim$10^-17 s^-1 and the AFGL 2591 value of $\sim$ $5.6 \times 10^{-17}$ s^-1; (2) depletion of CO at low temperatures ($\sim$20-40 K); (3) depletion of H₂CO and CH₃OH at moderate temperature (<60 K); and (4) variations in the initial abundances of a few other species. We discuss each of these differences below separately, as variations from the best fit model.

3.2 Effects of cosmic ray ionization rate

The ionization rate inferred in the modeling, $\sim$ $5 \times 10^{-16}$ s^-1, is much higher than the "standard'' cosmic ray ionization rate of $\sim$10^-17 s^-1 (e.g., Roberts & Herbst 2002; Black & Dalgarno 1977; O'Donnell & Watson 1974). To see the dependence of the model results on the ionization rate, Fig. 4 shows the mean difference $\Delta$ between the models and observations upon varying the ionization rate. The minimum deviation occurs in the range $\zeta = 5 \times 10^{-16}{-}10^{-15}$ s^-1. While the two values near the minimum are essentially indistinguishable, the ionization rate required for this fit is 50-100 times higher than the traditional cosmic-ray ionization rate used in dark cloud models (Lepp 1992).

Figure 4: Dependence of quality of fit as measured by the mean difference between predicted and observed line strengths, as a function of the ionization rate, $\zeta$. Notice the best fit near  $\zeta = 5\times 10^{-16}$ s-1. The two points near the minimum are indistinguishable at the level of uncertainty of the observations, and given the constraints of the model. In either case, the ionization rate is much larger than standard.

Open with DEXTER

The best fit range for $\zeta$ is expected to be meaningful, due to the fact that the mean difference is minimized here. The variation in $\Delta$ is damped by the fact that it is an average across all species. As a result, a $25\%$ variation in $\Delta$can correspond to a factor of 2 change in 1/4 of the species, or a factor of $\sim$8 change in $\sim$4-5 transitions. This is seen for the case of CO in Figs. 7 and 8 where a smaller change in $\Delta$ over all species corresponds to variation in physical parameters by $\sim$$100\times$, and variations in line strengths by up to $8 \times$. As a result, we infer that a minimum in $\Delta$ and a variation of $\sim$$20\%$is sufficient to draw conclusions, which implies that the preferred value of  $\zeta = 5 \times 10^{-16}{-}10^{-15}$ s^-1 is meaningfully different from other values tested.

The ionization rate $\zeta$ is assumed to be uniform throughout the source. Species that are most affected by the variation in the ionization rate are HCO⁺, HCN, SO, and H₂CO and show improvements of up to 100%. While the HCO⁺ abundance should be directly related to the ionization rate throughout the envelope, Doty et al. (2002) show that the remainder are predominantly active above 100 K. This implies that the ionization rate may be position dependent, with the most affected species in the interior.

The origin of this enhanced ionization is of physical interest. Recent models and measurements infer up to 2 orders of magnitude variation in the cosmic ray ionization rate (McCall et al. 2003; Liszt 2003; Doty et al. 2002; van der Tak & van Dishoeck 2000). It is difficult to understand such an extreme variation from source to source, implying that some other physical mechanism may produce or contribute to the ionization. This is especially true as the recent "high'' inferred values for $\zeta$ are for diffuse clouds, while dense cloud models have historically required much smaller cosmic-ray ionization rates near $1{-}3 \times 10^{-17}$ s^-1. In the case of IRAS 16293, we suggest the possibility that the enhanced ionization is due to X-rays produced by magnetic activity associated with accretion onto the protostars. This could both produce the exceptionally high inferred ionization rate, and preferentially affect the warmer interior species. The effects of a central X-ray source will be presented in a forthcoming paper (Doty et al., in preparation).

3.3 Effects of CO desorption temperature

	Figure 5: Dependence of quality of fit as measured by the mean difference between predicted and observed line strengths for the full chemistry / observational set, as a function of the CO desorption temperature ( $T_{{\rm CO}}$). Here, $x({\rm CO}) = 10^{-4}$. Notice the best fit for $T_{\rm CO} \sim 20$ K.
Open with DEXTER

Studies of solid CO on icy dust grains (Collings et al. 2003; Fraser et al. 2003; Galloway & Herbst 1994; Sandford & Allamandola 1993a; Nair & Adamson 1970) suggest that the bulk of the CO evaporates in a step-wise fashion between 20 and 70 K - depending upon whether it is trapped inside or lies on top of the ice - much lower than the temperature at which water ice desorbs from grains (Fraser et al. 2001). This is consistent with our best fit model. In order to test this, we have varied the CO desorption temperature from our baseline best-fit model. The results are shown in Fig. 5.

Clearly, the best fit requires a CO desorption temperature near $\sim$20 K and <60 K. A lower temperature both yields a worse fit to the observational data by overproducing the ¹³CO and C¹⁸O emission by 25% and 68% respectively at $T_{{\rm CO}}=10$ K, and is inconsistent with laboratory results. Much higher temperatures yield significantly worse fits to the data, underproducing both ¹³CO and C¹⁸O line fluxes by a factor of 50 by  $T_{{\rm CO}}=100$ K. The species most affected (aside from CO itself) are the cyanogens CN, HCN and HNC, and the CO ion-molecule byproducts HCO⁺, CS, and H₂CO, all of which show variations between 40-300%.

Physically, a low desorption temperature near 20 K would be an indication that a significant fraction of the CO is not intermixed with the H₂O in the grain mantle. A number of suggestions for differentiation in the ice have been made, including differentiation in the gas prior to adsorption, differentiation in the ice due to chemical and physical processing, and differentiated freeze-out during the cooling time behind a shock which has liberated the grain mantles (e.g., Schutte 1997; Bergin et al. 1999). While it is difficult to comment on the first two scenarios, it is doubtful that shock processing is the main cause for IRAS 16293-2422. In particular, the products of shock chemistry do not dominate the bulk of the envelope, and the small linewidths observed for many species further suggest that a large fraction of the volume of the gas is not shocked. Note that this analysis does not exclude that some fraction of the CO also evaporates at higher temperatures. In fact, Jørgensen et al. (2002) conclude from their analysis of the CO  $3 \rightarrow 2$ and  $2 \rightarrow 1$ isotopic lines in a sample of class 0 objects that some CO must still be frozen out at temperatures above 25 K.

	Figure 6: Dependence of quality of fit as measured by the mean difference between predicted and observed line strengths as a function of the CO fractional abundance for $T>T_{{\rm CO}}=20$ K. Notice the best fit for $x({\rm CO}) \sim 10^{-4}$.
Open with DEXTER

3.4 Effects of CO abundance

Previous authors (Schöier et al. 2002; Ceccarelli et al. 2000a) have inferred (constant) CO abundances in the range of  10^-5-10^-4. As a result the inner ( $T>T_{{\rm des,CO}}$) CO abundance is varied, keeping all other parameters the same as in our best fit model. The results are presented in Fig. 6.

Figure 7: Contours of the percentage difference between model and observed CO line strengths, for various values of CO desorption temperature ( $T_{{\rm CO}}$), abundance of cold CO ( $T<T_{{\rm CO}}$), and abundance of warm CO ( $T>T_{{\rm CO}}$) denoted by the different line types. The solid, dotted, and dashed lines correspond to  $x({\rm CO,warm})=10^{-4}$, $5\times 10^{-4}$, and $2\times 10^{-5}$ respectively. The top panel gives the results for C17O, and the bottom panel the average results for 13CO, C17O, and C18O. Notice that the C17O results strongly favor  $x({\rm CO,warm})=10^{-4}$, as do the results averaged over the CO isotopomers.

Open with DEXTER

As can be seen, CO abundances of  $1{-}4 \times 10^{-4}$ are preferred. Lower abundances produce too little C¹⁸O emission to be consistent with observations. The best fit - based upon CO data alone - is $x({\rm CO}) = 10^{-4}$where models match C¹⁸O observations to within 2%. By  $x({\rm CO}) = 5 \times 10^{-5}$ and  $4 \times 10^{-4}$, the discrepancy for the C¹⁸O lines reaches 37% and 95% respectively. The effect on $\Delta$ is smaller for two reasons: first other species are included in the mean difference, and second the comparison for HCO⁺ and CN are both improved as $x({\rm CO})$ increases.

The HCO⁺ and CN abundances both rely upon CO via ion-molecule reactions. Since this relation is much more indirect, the roles of HCO⁺ and CN abundances in fixing the CO abundance are discounted, and  $x({\rm CO}) \sim 10^{-4}$ is preferred. The constant abundance of $\sim$ $4\times 10^{-5}$ inferred by Schöier et al. (2002) in the empirical model can now be understood as a weighted average of a very low CO abundance at T<20 K and a higher abundance of $\sim$10^-4 at T>20 K.

We have also considerd the combined effects of CO desorption temperature, cold ( $T<T_{{\rm CO}}$) CO abundance, and warm ( $T>T_{{\rm CO}}$) CO abundance. The results are presented in Fig. 7 where we plot the mean percentage difference between model and observed line strengths for C¹⁷O (top panel), and all CO isotopomers (bottom panel). Warm CO abundances are denoted by the different line types, corresponding to  $x({\rm CO,warm}) = 10^{-4}, 5\times10^{-4}, 2\times10^{-5}$ respectively. As can be seen, the CO data strongly prefer $x({\rm CO,warm}) \sim 10^{-4}$, consistent with our results above. This is confirmed by the results in Fig. 8 where we plot contours of the mean difference, $\Delta$, over all observed lines. Again, models with $x({\rm CO,warm}) \sim 10^{-4}$ are preferred.

Figure 8: Contours of quality of fit as measured by the mean difference between predicted and observed line strengths over all species as a function of the CO desorption temperature ( $T_{{\rm CO}}$), abundance of cold CO ( $T<T_{{\rm CO}}$), and abundance of warm CO ( $T>T_{{\rm CO}}$) denoted by the different line types. The solid, dotted, and dashed lines correspond to $x({\rm CO,warm})=10^{-4}$, $5\times 10^{-4}$, and $2\times 10^{-5}$ respectively. Notice the agreement with Fig. 7 in constraining $16<T_{{\rm CO}}<30$ K, $x({\rm CO,cold})<10^{-5}$, and $x({\rm CO,warm}) \sim 10^{-4}$.

Open with DEXTER

Perhaps even more interestingly, the results for both the CO and general chemical network allow us to simultaneously constrain the desorption temperature and cold CO abundance. Taking uncertainties in the observational data of $\pm$30% suggests  $18 < T_{{\rm CO}} {\rm (K)} < 23$ based upon the 25% contour level in Fig. 7. Even if the uncertainties are larger, the results in Figs. 7 and 8 provide outer limits of $16 < T_{{\rm CO}} {\rm (K)} < 30$. There may be a potential region of degeneracy in Figs. 7 and 8 as the cold CO abundance appears to increase with increasing  $T_{{\rm CO}}$. It should be noted that the cold CO abundance is still $\sim$10^-5 in this case, a result that may be explained by significant evaporation near  $T_{{\rm CO}}$ and some partial/gradual evaporation of CO presumably in an H₂O matrix at higher temperatures. These results are consistent with both the cut along  $T_{{\rm CO}}$ and the laboratory results discussed above.

Finally, there does seem to be evidence of depletion in the cold gas for  $T<T_{{\rm CO}}$. The comparison for both the CO and overall set of observed species suggest a relatively firm upper limit of  $x({\rm CO,cold})<10^{-5}$. The regions of best fit appear to encompass values of 3-30 times less ( $3 \times 10^{-7}$ to $3 \times 10^{-6}$). The upper value of 10^-5 signifies a depletion of 90%, while the lower values correspond to 97% and 99% depletion respectively. While the exact level of depletion is uncertain, these results do confirm a significant sink of gas-phase CO - presumably as ices onto dust grains in the cool exterior. Such high levels of CO depletion are consistent with those found in cold pre-stellar cores (e.g., Bacmann et al. 2003) and the large abundances of deuterated molecules detected in the outer envelope of IRAS 16293-2422 (van Dishoeck et al. 1995; Loinard et al. 2000; Parise et al. 2002)

3.5 Effects of H₂CO and CH₃OH depletion temperature

	Figure 9: Dependence of quality of fit as measured by the mean difference between predicted and observed line strengths as a function of the temperature of H2CO and CH3OH desorption (TX). Notice the best fit for TX > 30 K.
Open with DEXTER

We also consider the effects of modifying the H₂CO and CH₃OH desorption temperatures, T_X. In the case of H₂CO which may be produced at some level in the gas phase (Doty et al. 2002), this corresponds to the temperature at which significant production occurs. The results are shown in Fig. 9, where the mean difference is plotted as a function of T_X.

To within the uncertainties of the observations and sophistication of the models, the results for T_X=30-100 K are indistinguishable. Consequently, the results in Fig. 9 suggest  T_X > 30 K. In the case of CH₃OH, which is almost certainly formed on the grain surface (e.g., Tielens & Hagen 1982; Blake et al. 1987; Doty et al. 2002), this is probably an indication of desorption. Desorption in this temperature range is consistent with current chemical understanding of solid ices. Sandford & Allamandola (1993b) measure a pure CH₃OH desorption temperature (in space) of 70-80 K. This is in keeping with the fact that CH₃OH should have a lower desorption temperature than water due to its weaker hydrogen bonding. For comparison, a Clausius-Clapeyron calculation which reproduces the water evaporation temperature well suggests an evaporation temperature for CH₃OH of $\sim$87 K (Alsindi et al. 2003). Likewise, these results are consistent with the empirical modeling of Schöier et al. (2002), who found $T_X({\rm CH}_{3}{\rm OH}) > 50$ K. As a result, it is encouraging that laboratory work, theoretical calculations, empirical models, and the results of this work are all in agreement with a CH₃OH desorption temperature of  $60 < T({\rm K}) < 100$.

In the case of H₂CO, the meaning of this temperature is less clear. The deviation between the observed and predicted H₂CO line strengths caused by the overproduction of H₂CO in the model do not significantly change as T_X is increased - from 49% at  $T_{{\rm des}}=10$ K to 44% at  $T_{{\rm des}}=100$ K. While Doty et al. (2002) suggested it was possible that gas-phase reactions could play a role in the H₂CO chemistry, they did not identify any that would cause a significant "jump'' in abundances in this temperature range. While it is difficult to directly constrain the H₂CO abundance, it is clear that the modeling here is insensitive to the amount of cold H₂CO, and that while no jump is required, a jump due to desorption is not ruled out - consistent with  T_X > 40 K as found by Schöier et al. (2002).

3.6 Outliers

As can be seen in Fig. 3, the three species OCS, C₂H, and HCS⁺ are outliers in our models. These species yield line strengths which diverge from the observations by factors of $\sim$30, 30, and 100 respectively. Some deviation due to radiative transfer, geometrical, and line of sight effects is to be expected. However, the differences for these three species are discrepant from the other species considered, implying that while our adopted physical/chemical structure is reasonable, our knowledge of the chemistry is lacking.

In the case of OCS, the chemistry is uncertain, and many of the reaction rates are estimates without significant laboratory study. As a result, it is not suprising that a discrepancy exists. For C₂H, UV photodissociation in the outer region - while included - is not significant as most of the C₂H is produced above 100 K. The dominant production is via recombination of C₂H₅⁺, and reaction of C₃H₂⁺ with O. Destruction is mainly through the neutral-neutral reaction with O. While a lower oxygen abundance can increase the C₂H abundance, it also decreases the abundances of the other important oxygen-bearing species such as SO, SO₂, and CH₃OH to below the observations. On the other hand, it is interesting to note that the destruction reaction with O is assumed to be temperature independent (MFW). The existence of a reaction barrier or a temperature dependent rate of collisions would both tend to decrease the destruction, which would have the effect of raising the C₂H abundance closer to the observed levels.

The third outlier, HCS⁺ follows an ionization balance with CS via dissociative recombination, and reactions with HCO⁺, H₃⁺, and H₃O⁺. However, raising $\zeta$ to  $5 \times 10^{-15}$ s^-1 only changes the HCS⁺ discrepancy by $\sim$3%. Such high values of the cosmic ray rate raise the most sensitive ion - HCO⁺ - to levels far above that observed. While the major destruction mechanism of dissociative recombination has a somewhat strong (T^-0.75) temperature dependence, this rate has been measured in the laboratory, and is considered to be accurate to within 25% by MFW. This leaves the production reactions of (HCO⁺, H₃⁺, and H₃O⁺) + CS as potential sources of uncertainty. Each of these rates are estimated. As such, it would be useful to measure them in the laboratory to confirm the rates adopted by MFW.

Finally, while HNC only differs by a factor of 5 and is thus not a major outlier, the production of HNC is still not fully understood (e.g. Rodgers & Charnley 2001; Charnley et al. 2001; Liszt & Lucas 2001). In our model, HNC is produced primarily through the dissociative recombination of HCNH⁺ and H₂NC⁺. In the recombination of HCNH⁺, where the branching fractions are taken to be 50% for CN, 25% for HCN, and 25% for HNC, the HNCH⁺ abundance is determined mostly by an ionization equilibrium in which the primary production paths are proton transfer between HNC and (HCO⁺, H₃O⁺, H₃⁺). These are also the dominant destruction paths for HNC. On the other hand, H₂NC⁺ is formed by C $^{+} + {\rm NH}_{3} \rightarrow {\rm H}_{2}{\rm NC}^{+}$, and dissociatively recombines to form HNC and CN in a 10:1 ratio. We note that while too much HNC is produced in our model, the CN abundance is somewhat low. This combination suggests that the adopted branching fractions for the dissociative recombination should perhaps be reinvestigated. Adopting the results of Talbi & Herbst (1998) for C $^{+} + {\rm NH}_{3} \rightarrow {\rm HCNH}^{+}$ has little effect, suggesting further concentration on H₂NC⁺Alternatively, some HNC destruction route may be missing in the networks.

3.7 Comparison with empirical model

	Figure 10: Radial abundance profiles for the species considered in the text for  $t=3 \times 10^{3}$ and $3 \times 10^{4}$ years. The results of the empirical model of Schöier et al. (2002) are given by the lines with arrows, for comparison.
Open with DEXTER

In principle, the empirical modeling approach should - with sufficient parameter variation - mimic the results of the full chemical modeling. In practice, it is difficult to parametrically vary all species in a sufficiently meaningful yet complete manner. As discussed in Sect. 1, the approach adopted by many authors is to treat abundances as either constant, or as piecewise constant with "jumps'' at appropriate temperatures. This was the approach taken in Schöier et al. (2002). It is interesting to compare the inferred abundance distributions from the empirical modeling from the more detailed results of the full chemical modeling.

Quantitatively, the Schöier et al. (2002) empirical abundances reproduce the observations with  $\Delta=0.24$, about half of our best fit model, $\Delta=0.51$. This is, however, not a suprise. In the Schöier et al. (2002) model, there are many more free parameters as the abundances for each of the observed species are varied both above and below the assumed desorption temperature. Furthermore, these variations are done without respect to constraints on the chemical network or evolution time. On the other hand, we directly specify the initial abundances for only three of the observed species (CO, H₂CO, and CH₃OH), and are constrained by the chemical network and its evolution. Furthermore, as discussed previously, the majority of the inputs to the chemical network are taken directly from agreement reached on a high-mass hot-core source, AFGL 2591. Taken together, these results are strongly encouraging as the chemical network comparison gives a good fit for significantly fewer direct parameter variations, confirms the proposed age of the source from physical evolution models, and directly tests the validity and extensibility of the chemical networks.

As a more direct comparion, the radial abundance profiles for the species discussed in this paper are presented in Fig. 10, where the two limiting times of $3 \times 10^{3}$ years and  $3 \times 10^{4}$ years are plotted. In general, the abundances inferred from the empirical modeling are grossly consistent with those from the more detailed full chemical modeling. For most species with constant abundances, the inferred abundances in mostly the outer envelope are equivalent to within a factor of 3-10 at our best-fit time (here taken to be that intermediate between the two limits). In some cases, e.g. CN, CS, and CH₂CO, the chemical model abundance oscillates with position in the cloud and the empirical value is simply a rough average of these complicated profiles. The most significant discrepancies are for OCS, HCS⁺, C₂H, and HNC, which are discussed above.

Even more interesting is the comparison of the jump models. In general, those species which show significant spatial variation in the full chemical model are represented as "jumps'' in the Schöier et al. (2002) empirical model. Of these, the empirical and full chemical models generally agree to within a factor of 3 or so. As discussed in Sect. 3.4, the inferred CO abundance can be understood in such a "jump'' model. The significant discrepancy is H₂CO, which the full chemical model predicts to have only a small jump, while the empirical model infers a cold, outer abundance some 15 times lower.

3.8 Predictions for future observations

	Figure 11: Radial abundance profiles for the some of the more abundant species that may be targets for future observations. The data are reproduced for $t=3 \times 10^{3}$ and $3 \times 10^{4}$ years.
Open with DEXTER

While the model presented is able to simultaneously match the SED and many of the molecular line observations, a significant test will be the predictions it makes that can be studied by future facilities. In particular, with CARMA, SIRTF, and SOFIA data to be available in the next few years, and the upcoming leaps in resolution and sensitivity from Herschel and ALMA, it will be possible to probe many of the transitions and much of the spatial structure of IRAS 16293.

To aid such future observations, the radial distributions of a number of interesting species are shown in Fig. 11. Also, column densities predicted by the model toward IRAS 16293-2422 for $t=3 \times 10^{3}$ years and  $3 \times 10^{4}$ years are given in Table 4. This includes, in particular, radial column densities given by $N_{X,{\rm rad}} = \int n_{X}(r) {\rm d}r$, where n_X(r) is the density of species X as a function of r. We also include the column density averaged over a 15 arcsec beam toward IRAS 16293, given by  $N_{X,{\rm beam}} = \int \int n_{X}(z,b) {\rm d}z G(b) 2\pi b {\rm d}b / \int G(b) 2 \pi b {\rm d}b$, where b is the impact parameter, and G(b) is the beam response function. The column densities are sorted from highest to lowest in the 15 arcsec beam at  $t=3 \times 10^{3}$ years, and continue down to $\sim$10¹³ cm^-2.

Table 4: IRAS 16293-2422 predicted column densities (cm^-2).

Finally, it is intruiging to speculate on the effects and potential observability of a collapsing cloud. The evolution of a gas parcel in a protostellar envelope passing from a low temperature/density to a high temperature/density through infall has been studied by CHT96, and Rodgers & Charnley (2003) among others. In the outer regions, cool ion-molecule chemistry dominates. As the gas heats while infalling, each adsorbed species passes through a sublimation front. In the interior warm neutral chemistry can take place, for example the conversion of gas-phase oxygen to H₂O (e.g., CHT96; Charnley 1997), and the production of significant complex cyanogens and hydrocarbons (e.g., Rodgers & Charnley 2001; Doty et al. 2002). In free-fall collapse of low-mass YSOs, the dynamical timescale in the warm interior is less than the chemical timescale, leading to inner core abundances that mirror those in the cool exterior. The observations of a wide variety of complex daughter species observed in these warm interiors (e.g. Cazaux et al. 2003) can only be understood if the collapse/infall is "slowed'' to timescales over at least 10⁴ years (Rodgers & Charnley 2003) - a scenario more in-line with our adopted static case. Since the chemistry encodes the temperature/density temporal evolution of the gas, small-scale spatial variations in abundances should, in principle, be able to distinguish between various dynamical scenarios such as static, Shu (1977) collapse, Larson-Penston (Larson 1969; Penston 1969) infall, etc. This will, however, require a next generation of modeling that includes detailed dynamics, thermal balance, chemistry, and radiative transfer, as well as high spatial and spectral resolution observations with instruments such as Herschel, ALMA, CARMA, and SOFIA to probe multiple lines at $\sim$100 AU resolution.

4 Conclusions

We have constructed detailed thermal and gas-phase chemical models for IRAS 16293-2422 based upon the physical model of Schöier et al. (2002). These models were used to probe the validity of the proposed physical structure, as well as study the chemical evolution of the source, and to test the application of our combined "hot-core''/envelope chemistry model of AFGL 2591 to a low-mass "hot-core''-like source. In particular, we find that:

1.

The combined application of a physical, thermal, and chemical model with detailed radiative transfer is a powerful tool in constraining the structure and evolution of depth-dependent sources.

2.

The time and position dependent model of Doty et al. (2002) can be meaningfully applied to IRAS 16293-2422, yielding results qualitatively similar to the massive YSO AFGL 2591.

3.

The best fit for IRAS 16293-2422 occurs for times in the range $3 \times 10^{3} < t({\rm yrs}) < 3 \times 10^{4}$, consistent with existing infall models (Sect. 3.1).

4.

The best-fit ionization rate in IRAS 16293-2422 is high, $\sim$ $5 \times 10^{-16}{-}10^{-15}$ s^-1, compared with previous results for dense clouds. We propose that this may be due to X-ray emission from the central sources (Sect. 3.2).

5.

Our best fit suggests that important CO desorption occurs at low temperatures, $\sim$20 K, and constrained to the range $15 < T_{{\rm CO}} {\rm (K)} < 40$. Some solid CO may remain at higher temperatures for these timescales. These results are in agreement with recent laboratory data of CO on - but not mixed with - a water ice (Sect. 3.3).

6.

We can also constrain the warm ( $T>T_{{\rm CO}}$) and cold ( $T<T_{{\rm CO}}$) CO abundances to $x({\rm CO,warm}) \sim 10^{-4}$ and  $x({\rm CO,cold})<10^{-5}$ and most probably <10^-6. These results are reflected in both the CO lines and in the results from the greater chemical network, and suggest significant (>90%) depletions at low temperatures (Sect. 3.4).

7.

CH₃OH appears to desorb at temperatures $60 < T({\rm K}) < 100$, consistent with laboratory expectations. On the other hand, the comparison between the H₂CO predictions and observations are insensitive to the amount of cold H₂CO present (Sect. 3.5).

8.

The chemistry of HNC, C₂H, and HCS⁺ may not be fully understood. In particular, it would be useful to measure the branching fraction of dissociative recombination of H₂NC⁺, the temperature dependence of the reaction O+C₂H, and the ion-molecule reaction rates (HCO⁺, H₃⁺, H₃O⁺) + CS (Sect. 3.6).

Acknowledgements

The authors are grateful to Jes Jørgensen and Helen Fraser for fruitful discussions. We thank the referee for comments which improved the manuscript. This work was partially supported under grants from The Research Corporation (SDD), and the Netherlands Organisation for Scientific Research (NWO) through grant 614.041.004. F.L.S. further acknowledges financial support from the Swedish Reseach Council. Astrochemistry at Leiden is supported through an NWO Spinoza award.

References

• Alsindi, W. Z., McCoustra, M. R. S., George, M. W., & Fraser, H. J. 2003, in preparation (In the text)
• Bacmann, A., Lefloch, B., Ceccarelli, C., et al. 2003, ApJ, 585, L55 [NASA ADS] [CrossRef] (In the text)
• Bergin, E. A., Neufeld, D. A., & Melnick, G. J. 1999, ApJ, 510, L145 [NASA ADS] [CrossRef] (In the text)
• Black, J. H., & Dalgarno, A. 1977, ApJS, 34, 405 [NASA ADS] [CrossRef] (In the text)
• Blake, G. A., Sutton, E. C., Masson, C. R., & Phillips, T. G. 1987, ApJ, 315, 621 [NASA ADS] [CrossRef] (In the text)
• Blake, G. A., van Dishoeck, E. F., Jansen, D. J., Groesbeck, T. D., & Mundy, L. G. 1994, ApJ, 428, 680 [NASA ADS] [CrossRef] (In the text)
• Boonman, A. M. S., van Dishoeck, E. F., Lahuis, F., Wright, C. M., & Doty, S. D. 2000, in ISO beyond the peaks, ed. A. Salama, M. F. Kessler, K. Leech, & B. Schulz, ESA-SP 456, 67 [arXiv:astro-ph/0105249]
• Boonman, A. M. S., Doty, S. D., van Dishoeck, E. F., et al. 2003, A&A, 406, 937 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Brown, P. D., & Charnley, S. B. 1990, MNRAS, 244, 432 [NASA ADS] (In the text)
• Castets, A., Ceccarelli, C., Loinard, L., Caux, E., & Lefloch, B. 2001, A&A, 375, 40 [EDP Sciences] [NASA ADS] [CrossRef]
• Cazaux, S., Tielens, A. G. G. M., Ceccarelli, C., et al. 2003, ApJ, 593, L51 [NASA ADS] [CrossRef] (In the text)
• Ceccarelli, C., Castets, A., Caux, E., et al. 2000a, A&A, 355, 1129 [NASA ADS] (In the text)
• Ceccarelli, C., Hollenbach, D. J., & Tielens, A. G. G. M. 1996, ApJ, 471, 400 [NASA ADS] [CrossRef] (In the text)
• Ceccarelli, C., Loinard, L., Castets, A., Tielens, A. G. G. M., & Caux, E. 2000b, A&A, 357, L9 [NASA ADS] (In the text)
• Charnley, S. B. 1997, ApJ, 481, 396 [NASA ADS] [CrossRef] (In the text)
• Charnley, S. B., Rodgers, S. D., & Ehrenfruend, P. 2001, A&A, 378, 1024 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Collings, M. P., Dever, J. W., Fraser, H. J., McCoustra, M. R. S., & Wiliams, D. A. 2003, ApJ, 583, 1058 [NASA ADS] [CrossRef] (In the text)
• Doty, S. D., & Neufeld, D. A. 1997, ApJ, 489, 122 [NASA ADS] [CrossRef] (In the text)
• Doty, S. D., & Palotti, M. L. 2002, MNRAS, 335, 993 [NASA ADS] [CrossRef] (In the text)
• Doty, S. D., Stäuber, P., van Dishoeck, E. F., & Schöier, F. L. 2003, A&A, in preparation
• Doty, S. D., van Dishoeck, E. F., van der Tak, F. F. S., & Boonman, A. M. S. 2002, A&A, 389, 446 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Evans, N. J. 1999, ARA&A, 37, 311 [NASA ADS] (In the text)
• Fraser, H. J., Collings, M. P., McCoustra, M. R. S., & Williams, D. A. 2001, MNRAS, 327, 1165 [NASA ADS] [CrossRef] (In the text)
• Fraser, H. J., Collings, M. P., McCoustra, M. R. S., & Dever, J. W. 2003, in preparation (In the text)
• Galloway, E. T., & Herbst, E. 1994, A&A, 287, 633 [NASA ADS] (In the text)
• Ivezic, Z., Nenkova, M., & Elitzur, M. 1999, User Manual for DUSTY (Univ. Kentucky Internal Rep.) (In the text)
• Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2002, A&A, 389, 908 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Larson, R. B. 1969, MNRAS, 145, 271 [NASA ADS] (In the text)
• Lepp, S. 1992, in Astrochemistry of Cosmic Phenomena, ed. P. D. Singh (Dordrecht: Kluwer), Proc. IAU Symp., 150, 471 (In the text)
• Liszt, H. 2003, A&A, 398, 621 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Liszt, H., & Lucas, R. 2001, A&A, 370, 576 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Loinard, L., Castets, A., Ceccarelli, C., et al. 2000, A&A, 359, 1169 [NASA ADS] (In the text)
• Looney, L. W., Mundy, L. G., & Welch, W. J. 2000, ApJ, 529, 477 [NASA ADS] [CrossRef] (In the text)
• McCall, B. J., Huneycutt, A. J., Saykally, R. J., et al. 2003, Nature, 422, 500 [NASA ADS] [CrossRef] (In the text)
• Millar, T., Farquhar, P., & Willacy, K. (MFW) 1997, A&AS, 121, 193 (In the text)
• Millar, T. J., & Freeman, A. 1984, MNRAS, 207, 405 [NASA ADS] (In the text)
• Millar, T. J., MacDonald, G. H., & Gibb, A. G. 1997, A&A, 325, 1163 [NASA ADS] (In the text)
• Millar, T. J., & Nejad, L. A. M. 1985, MNRAS, 217, 507 [NASA ADS]
• Mundy, L. G., Wootten, A., Wilking, B. A., Blake, G. A., & Sargent, A. I. 1992, ApJ, 385, 306 [NASA ADS] [CrossRef]
• Nair, N. K., & Adamson, A. W. 1970, J. Phys. Chem., 74, 2229 [NASA ADS] (In the text)
• O'Donnell, E. F., & Watson, W. D. 1974, ApJ, 191, 89 [NASA ADS] [CrossRef] (In the text)
• Parise, B., Ceccarelli, C., Tielens, A. G. G. M., et al. 2002, A&A, 393, L49 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Penston, M. V. 1969, MNRAS, 144, 425 [NASA ADS] (In the text)
• Roberts, H., & Herbst, E. 2002, A&A, 395, 233 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Rodgers, S. D., & Charnley, S. B. 2001, ApJ, 546, 324 [NASA ADS] [CrossRef] (In the text)
• Rodgers, S. D., & Charnley, S. B. 2003, ApJ, 585, 355 [NASA ADS] [CrossRef] (In the text)
• Sandford, S. A., & Allamandola, L. J. 1993a, Icarus, 106, 478 [NASA ADS] [CrossRef] (In the text)
• Sandford, S. A., & Allamandola, L. J. 1993b, ApJ, 417, 815 [NASA ADS] [CrossRef] (In the text)
• Schöier, F. L. 2000, Ph.D. Thesis, Stockholm Observatory, Sweden (In the text)
• Schöier, F. L., Jørgensen, J. K., van Dishoeck, E. F., & Blake, G. A. 2003, submitted (In the text)
• Schöier, F. L., Jørgensen, J. K., van Dishoeck, E. F., & Blake, G. A. 2002, A&A, 390, 1001 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Schutte, W. A. 1997, in The Cosmic Dust Connection, Proc. of the NATO Advanced Study Institute, ed. J. M. Greenberg (Dordrecht: Reidel), 1 (In the text)
• Shirley, Y. L., Evan, N. J., Rawlings, J. M. C., & Gregersen, E. M. 2000, ApJS, 131, 249 [NASA ADS] [CrossRef] (In the text)
• Shu, F. H. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] (In the text)
• Shu, F., Najita, J., Galli, D., Ostriker, E., & Lizano, S. 1993, in Protostars and Planets III, ed. E. H. Levy, & J. Lunine (Tucson: Univ. Arizona Press), 3 (In the text)
• Talbi, D., & Herbst, E. 1998, A&A, 333, 1007 [NASA ADS] (In the text)
• Terzieva, R., & Herbst, E. 1998, ApJ, 501, 270 [NASA ADS] [CrossRef] (In the text)
• Tielens, A. G. G. M., & Hagen, W. 1982, A&A, 114, 245 [NASA ADS] (In the text)
• van der Tak, F. F. S., & van Dishoeck, E. F. 2000, A&A, 358, L79 [NASA ADS] (In the text)
• van der Tak, F. F. S., van Dishoeck, E. F., Evans, N. J., Bakker, E. J., & Blake, G. A. 1999, ApJ, 522, 991 [NASA ADS] [CrossRef]
• van der Tak, F. F. S., van Dishoeck, E. F., Evans, N. J., & Blake, G. A. 2000, ApJ, 537, 283 [NASA ADS] [CrossRef] (In the text)
• van Dishoeck, E. F. 2003, to appear in Chemistry as a Diagnostic of Star Formation, ed. C. L. Curry, & M. Fich (Ottawa: NRC Press) (In the text)
• van Dishoeck, E. F., Blake, G. A., Jansen, D. J., & Groesbeck, T. D. 1995, ApJ, 447, 760 [NASA ADS] [CrossRef] (In the text)
• van Dishoeck, E. F., & van der Tak, F. F. S. 2000, in Astrochemistry: From Molecular Clouds to Planetary Systems, ed. Y. C. Minh, & E. F. van Dishoeck (San Francisco: Astron. Soc. Pac.), Proc. IAU Symp., 197, 97 (In the text)
• van Zadelhoff, G.-J., Dullemon, C. P., van der Tak, F. F. S., et al. 2002, A&A, 395, 373 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Viti, S., Natarajan, S., & Williams, D. A. 2002, MNRAS, 336, 797 [NASA ADS] [CrossRef]
• Walmsley, C. M., & Schilke, P. 1992, in Astrochemistry of Cosmic Phenomena, ed. P. D. Singh (Dordrecht: Kluwer), Proc. IAU Symp., 150, 251 (In the text)
• Whittet, D. C. B. 1974, MNRAS, 168, 371 [NASA ADS] (In the text)