A&A 418, 1021-1034 (2004)
S. D. Doty 1 - F. L. Schöier 2,3 - E. F. van Dishoeck 3
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
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 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 90% or more of the CO is frozen out. The inability of the model to explain the HCS+, C2H, 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
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
as a function of radius r and the dust temperature
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''
is taken to be equal to
and the excitation, radiative transfer and fluxes of the various
molecular lines are calculated for an assumed abundance profile
(H2)(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 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 (160 pc, Whittet 1974) low mass, low luminosity (27 ) protostellar object located within the 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 H2CO, H2O, and SiO in a "full chemical model'', and found strong evidence for large abundance enhancements of these species in the innermost part (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.
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 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 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 H2CO 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 30%.
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.
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 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 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 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).
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., H2O). The only exceptions are CO, which we take to desorb at , and H2CO and CH3OH, which we take to desorb or "jump'' at TX.
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 H2 for AFGL 2591.
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 500 K in the ground vibrational state of the observed molecules. Chemical evolution times from years to years are considered, with one dex spacing.
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,
The parameters varied in the models are the cosmic ray ionization rate , the adopted initial abundances in the inner and outer regions, and the desorption temperatures of selected species (CO, H2CO, CH3OH). Detailed radial profiles of selected species are presented in Sects. 3.7 and 3.8.
|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 . 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 , 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 .
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 deviation in 13CO, and the small uncertainties on 13CO, C18O, and HCO+. It is possible that the 13CO discrepancy could be due to a different 12C/13C 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 .
In the comparison there exist three species which deviate by more than an order of magnitude. The outliers are: OCS, C2H, 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 when these are omitted. Including them raises 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 10-17 s-1 and the AFGL 2591 value of s-1; (2) depletion of CO at low temperatures (20-40 K); (3) depletion of H2CO and CH3OH 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.
The ionization rate inferred in the modeling, s-1, is much higher than the "standard'' cosmic ray ionization rate of 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 between the models and observations upon varying the ionization rate. The minimum deviation occurs in the range 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, . Notice the best fit near 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 is expected to be meaningful, due to the fact that the mean difference is minimized here. The variation in is damped by the fact that it is an average across all species. As a result, a variation in can correspond to a factor of 2 change in 1/4 of the species, or a factor of 8 change in 4-5 transitions. This is seen for the case of CO in Figs. 7 and 8 where a smaller change in over all species corresponds to variation in physical parameters by , and variations in line strengths by up to . As a result, we infer that a minimum in and a variation of is sufficient to draw conclusions, which implies that the preferred value of s-1 is meaningfully different from other values tested.
The ionization rate 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 H2CO 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 are for diffuse clouds, while dense cloud models have historically required much smaller cosmic-ray ionization rates near 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).
|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 ( ). Here, . Notice the best fit for 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 20 K and <60 K. A lower temperature both yields a worse fit to the observational data by overproducing the 13CO and C18O emission by 25% and 68% respectively at K, and is inconsistent with laboratory results. Much higher temperatures yield significantly worse fits to the data, underproducing both 13CO and C18O line fluxes by a factor of 50 by 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 H2CO, 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 H2O 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 and 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 K. Notice the best fit for .|
|Open with DEXTER|
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 ( ) 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 ( ), abundance of cold CO ( ), and abundance of warm CO ( ) denoted by the different line types. The solid, dotted, and dashed lines correspond to , , and 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 , as do the results averaged over the CO isotopomers.|
|Open with DEXTER|
As can be seen, CO abundances of are preferred. Lower abundances produce too little C18O emission to be consistent with observations. The best fit - based upon CO data alone - is where models match C18O observations to within 2%. By and , the discrepancy for the C18O lines reaches 37% and 95% respectively. The effect on 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 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 is preferred. The constant abundance of 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 10-4 at T>20 K.
We have also considerd the combined effects of CO desorption
temperature, cold (
CO abundance, and warm (
CO abundance. The results are presented in
Fig. 7 where we plot the mean percentage difference
between model and observed line strengths for C17O (top panel),
and all CO isotopomers (bottom panel). Warm CO abundances are denoted
by the different line types, corresponding to
respectively. As can be
seen, the CO data strongly prefer
consistent with our results above. This is confirmed
by the results in Fig. 8 where we plot contours of the
mean difference, ,
over all observed lines. Again, models
|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 ( ), abundance of cold CO ( ), and abundance of warm CO ( ) denoted by the different line types. The solid, dotted, and dashed lines correspond to , , and respectively. Notice the agreement with Fig. 7 in constraining K, , and .|
|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 30% suggests 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 . There may be a potential region of degeneracy in Figs. 7 and 8 as the cold CO abundance appears to increase with increasing . It should be noted that the cold CO abundance is still 10-5 in this case, a result that may be explained by significant evaporation near and some partial/gradual evaporation of CO presumably in an H2O matrix at higher temperatures. These results are consistent with both the cut along and the laboratory results discussed above.
Finally, there does seem to be evidence of depletion in the cold gas for . The comparison for both the CO and overall set of observed species suggest a relatively firm upper limit of . The regions of best fit appear to encompass values of 3-30 times less ( to ). 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)
|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 H2CO and CH3OH desorption temperatures, TX. In the case of H2CO 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 TX.
To within the uncertainties of the observations and sophistication of the models, the results for TX=30-100 K are indistinguishable. Consequently, the results in Fig. 9 suggest TX > 30 K. In the case of CH3OH, 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 CH3OH desorption temperature (in space) of 70-80 K. This is in keeping with the fact that CH3OH 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 CH3OH of 87 K (Alsindi et al. 2003). Likewise, these results are consistent with the empirical modeling of Schöier et al. (2002), who found 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 CH3OH desorption temperature of .
In the case of H2CO, the meaning of this temperature is less clear. The deviation between the observed and predicted H2CO line strengths caused by the overproduction of H2CO in the model do not significantly change as TX is increased - from 49% at K to 44% at K. While Doty et al. (2002) suggested it was possible that gas-phase reactions could play a role in the H2CO 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 H2CO abundance, it is clear that the modeling here is insensitive to the amount of cold H2CO, and that while no jump is required, a jump due to desorption is not ruled out - consistent with TX > 40 K as found by Schöier et al. (2002).
As can be seen in Fig. 3, the three species OCS, C2H, and HCS+ are outliers in our models. These species yield line strengths which diverge from the observations by factors of 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 C2H, UV photodissociation in the outer region - while included - is not significant as most of the C2H is produced above 100 K. The dominant production is via recombination of C2H5+, and reaction of C3H2+ with O. Destruction is mainly through the neutral-neutral reaction with O. While a lower oxygen abundance can increase the C2H abundance, it also decreases the abundances of the other important oxygen-bearing species such as SO, SO2, and CH3OH 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 C2H abundance closer to the observed levels.
The third outlier, HCS+ follows an ionization balance with CS via dissociative recombination, and reactions with HCO+, H3+, and H3O+. However, raising to s-1 only changes the HCS+ discrepancy by 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+, H3+, and H3O+) + 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 H2NC+. 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+, H3O+, H3+). These are also the dominant destruction paths for HNC. On the other hand, H2NC+ is formed by C , 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 has little effect, suggesting further concentration on H2NC+Alternatively, some HNC destruction route may be missing in the networks.
|Figure 10: Radial abundance profiles for the species considered in the text for and 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 , about half of our best fit model, . 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, H2CO, and CH3OH), 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 years and 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 CH2CO, 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+, C2H, 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 H2CO, 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.
|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 and 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 years and years are given in Table 4. This includes, in particular, radial column densities given by , where nX(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 , 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 years, and continue down to 1013 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 H2O (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 104 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 100 AU resolution.
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.