A&A 389, 908-930 (2002)
J. K. Jørgensen - F. L. Schöier - E. F. van Dishoeck
Leiden Observatory, PO Box 9513, 2300 RA Leiden, The Netherlands
Received 14 December 2001 / Accepted 30 April 2002
We present 1D radiative transfer modelling of the envelopes of a sample of 18 low-mass protostars and pre-stellar cores with the aim of setting up realistic physical models, for use in a chemical description of the sources. The density and temperature profiles of the envelopes are constrained from their radial profiles obtained from SCUBA maps at 450 and 850 m and from measurements of the source fluxes ranging from 60 m to 1.3 mm. The densities of the envelopes within 10000 AU can be described by single power-laws for the class 0 and I sources with ranging from 1.3 to 1.9, with typical uncertainties of 0.2. Four sources have flatter profiles, either due to asymmetries or to the presence of an outer constant density region. No significant difference is found between class 0 and I sources. The power-law fits fail for the pre-stellar cores, supporting recent results that such cores do not have a central source of heating. The derived physical models are used as input for Monte Carlo modelling of submillimeter C18O and C17O emission. It is found that class I objects typically show CO abundances close to those found in local molecular clouds, but that class 0 sources and pre-stellar cores show lower abundances by almost an order of magnitude implying that significant depletion occurs for the early phases of star formation. While the 2-1 and 3-2 isotopic lines can be fitted using a constant fractional CO abundance throughout the envelope, the 1-0 lines are significantly underestimated, possibly due to contribution of ambient molecular cloud material to the observed emission. The difference between the class 0 and I objects may be related to the properties of the CO ices.
Key words: stars: formation - ISM: molecules - ISM: abundances - stars: circumstellar matter - radiative transfer - astrochemistry
The Shu (1977) model predicts that the outer parts of the envelope follow a density profile similar to the solution of an isothermal sphere (Larson 1969), while within the collapse radius, which is determined by the sound speed in the envelope material and the time since the onset of the collapse, the density tends to flatten nearing a profile in the innermost parts. This model has subsequently been refined to include for example rotational flattening (Terebey et al. 1984) where such effects are described as a perturbation to the Shu (1977) solution.
An open question about the properties of YSOs in the earliest protostellar stage is how the structure of the envelope reflects the initial protostellar collapse and how it will affect the subsequent evolution of the protostar, for example in defining the properties of the circumstellar disk from which planets may be formed later. One of the possible shortcomings of the Shu-model is the adopted, constant, accretion rate. Foster & Chevalier (1993) performed hydrodynamical simulations of the stages before the protostellar collapse and found that the structure of the core at the collapse initiating point is highly dependent on the initial conditions; only in the case where a large ratio exists between the radius of an outer envelope with a flat density profile and an inner core with a steep density profile, will the core evolve to reproduce the conditions in the Shu model. In an analytical study, Henriksen et al. (1997) suggested that the accretion history of protostars could be divided into two phases for cores with a flat inner density profile: a violent early phase with high accretion rates (corresponding to the class 0 phase) that declines until a phase with mass accretion rates similar to the predictions in the Shu-model is reached (class I objects), i.e., a distinction between class 0 and I objects based on ages. Whether this is indeed the case has recently been questioned by Jayawardhana et al. (2001), who instead suggest that both class 0 and I objects are protostellar in nature, but just associated with environments of different physical properties, with the class 0 objects in more dense environments leading to the higher accretion rates observed towards these sources.
The chemical composition of the envelope may be an alternative tracer of the evolution. Indeed, for high-mass YSOs, combined infrared and submillimeter data have shown systematic heating trends reflected in the ice spectra, gas/ice ratios and gas-phase abundances (e.g., Gerakines et al. 1999; Boogert et al. 2000; van der Tak et al. 2000a). One of the prime motivations for this work is to extend similar chemical studies to low-mass objects, and extensive (sub-)millimeter line data for a sample of such sources are being collected at various telescopes, which can be complemented by future SIRTF infrared data.
In order to address these issues the physical parameters within the envelope, in particular the density and temperature profiles and the velocity field, are needed. The first two can be obtained through modelling of the dust continuum emission observed towards the sources, while observations of molecules like CO and CS can trace the gas component and velocity field. At the densities observed in the inner parts of the envelopes of YSOs it is reasonable to expect gas and dust coupling, which is usually expressed by the canonical dust-to-gas ratio of 1:100 and the assumption that the dust and gas temperatures are similar. Therefore a physical model for the envelopes derived on the basis of the dust emission can be used as input for modelling of the abundances of the various molecules.
|(hh mm ss)||(dd mm ss)||(K)||()||(pc)|
|L1448-I2||03 25 22.4||+30 45 12||60||3||220||03222+3034||Class 0|
|L1448-C||03 25 38.8||+30 44 05||54||5||220||L1448-MM|
|N1333-I2||03 28 55.4||+31 14 35||50||16||220||03258+3104/SVS19|
|N1333-I4A||03 29 10.3||+31 13 31||34||6||220|
|N1333-I4B||03 29 12.0||+31 13 09||36||6||220|
|L1527||04 39 53.9||+26 03 10||36||2||140||04368+2557|
|VLA1623||16 26 26.4||-24 24 30||35||1||160|
|L483||18 17 29.8||-04 39 38||52||9||200||18148-0440|
|L723||19 17 53.7||+19 12 20||47||3||300||19156+1906|
|L1157||20 39 06.2||+68 02 22||42||6||325||20386+6751|
|CB244||23 25 46.7||+74 17 37||56||1||180||23238+7401/L1262|
|IRAS 16293-2422||16 32 22.7||-24 28 32||43||27||160||(a)|
|L1489||04 04 43.0||+26 18 57||238||3.7||140||04016+2610||Class I|
|TMR1||04 39 13.7||+25 53 21||144||3.7||140||04361+2547|
|L1551-I5||04 31 34.1||+18 08 05||97||28||140||04287+1801||(b)|
|TMC1A||04 39 34.9||+25 41 45||172||2.2||140||04365+2535||(b)|
|TMC1||04 41 12.4||+25 46 36||139||0.66||140||04381+2540||(b)|
|L1544||05 04 17.2||+25 10 44||18||1||140||Pre-stellar|
|L1689B||16 34 49.1||-24 37 55||18||0.2||160|
To further explore these properties of the protostellar envelopes we have undertaken full 1D radiative transfer modelling of a sample of protostars and pre-stellar cores (see Sect. 2.1) using the radiative transfer code DUSTY (Ivezic & Elitzur 1997). Assuming power-law density distributions we solve for the temperature distribution and constrain the physical parameters of the envelopes by comparison of the results from the modelling to SCUBA images of the individual sources and their spectral energy distributions (SEDs) using a rigorous method. Besides giving a description of the physical properties of low-mass protostellar envelopes, the derived density and temperature profiles are essential as input for detailed chemical modelling of molecules observed towards these objects. Also, a good description of the envelope structure is needed to constrain the properties of the disks in the embedded phase (e.g., Keene & Masson 1990; Hogerheijde et al. 1998,1999; Looney et al. 2000).
In Sect. 2 our sample of sources is presented and the reduction and calibration briefly discussed (see also Schöier et al. 2002). In Sect. 3 the modelling of the sources is described and the derived envelope parameters presented. The properties of the individual sources are described in Sect. 3.4. In Sect. 4 the implication of these results are discussed and compared to other work done in this field. The results of the continuum modelling will be used in a later paper as physical input for detailed radiative transfer modelling of molecular line emission for the class 0 objects - as has been done for class I objects (Hogerheijde et al. 1998) and high-mass YSOs (van der Tak et al. 2000b) and as presented for the low-mass class 0 object IRAS 16293-2422 (Ceccarelli et al. 2000a,b; Schöier et al. 2002). In Sect. 5, the first results of this radiative transfer analysis for C18O and C17O is presented.
For the class 0 objects we have adopted luminosities and distances from André et al. (2000), for the class I objects the values from Motte & André (2001) and for the pre-stellar cores and CB244 distances and luminosities from Shirley et al. (2000). There are a few exceptions, however: for the objects related to the Perseus region we assume a distance of 220 pc and scale the luminosities from André et al. (2000) accordingly, while a distance of 325 pc is assumed for L1157 as in Shirley et al. (2000). The sample is summarized in Table 1. The class 0 object IRAS 16293-2422 treated in Schöier et al. (2002) has been included for comparison here as well.
For each source the flux scale was calibrated using available data for one of the standard calibrators, either a planet or a strong submillimeter source like CRL618. From the calibrated maps the total integrated fluxes were derived and the 1D brightness profiles were extracted by measuring the flux in annuli around the peak flux. The annuli were chosen with radii of half the beam (4.5 for the 450 m data and 7.5 for the 850 m data) so that a reasonable noise-level is obtained, while still making the annuli narrow enough to get information about the source structure without oversampling the data. Actually the spread in the fluxes measured for the points in each annulus due to instrumental and calibration noise was negligible compared to the spread due to (1) the gradient in brightness across each annuli, and (2) deviations from circular structure of the sources. One problem in extracting the brightness profiles was presented by cases where nearby companions were contributing significantly when complete circular annuli were constructed - the most extreme example being N1333-I4 with two close protostars. In these cases emission from "secondary'' components was blocked out by simply not including data-points in the direction of these closeby sources when calculating the mean flux in each annulus.
The derived parameters for L1157 and CB244 are given
in Table 3. Images of the two sources at the two SCUBA
wavelengths are presented in Fig. 1. As seen from the
figure, both sources are quite circular with only a small degree of
extended emission. Comparison with the 850 m data of
Shirley et al. (2000) for the 40
aperture shows that the fluxes
agree well within the 20% uncertainty assumed for the calibration.
|450 m||850 m||450 m||850 m|
|Figure 1: SCUBA images of CB244 and L1157 at 450 and 850 m. The contours indicate the intensity corresponding to 2, 4, etc. with being the RMS noise given for each source and wavelength in Table 3.|
The integrated line data were brought from the antenna temperature scale to the main-beam brightness scale by dividing by the main-beam brightness efficiency taken to be 0.69 for data obtained using the JCMT A band receivers (210-270 GHz; the J=2-1 transitions) and 0.59-0.63 for respectively the old B3i (before December 1996) and new B3 receivers (330-370 GHz; the J=3-2 transitions). For the IRAM 30 m observations beam efficiencies of 0.74 and 0.54 and forward efficiencies of 0.95 and 0.91 were adopted for respectively the C17O J=1-0 and J=2-1 lines, which corresponds to main-beam brightness efficiencies ( ) of respectively 0.78 and 0.59. For the Onsala 20 m telescope was adopted for the C18O and C17O J=1-0lines. The relevant beam sizes for the JCMT are 21 and 14 at respectively 220 and 330 GHz, for the IRAM 30 m, 22 and 11 at respectively 112 and 224 GHz and for the Onsala 20 m, 33. The velocity resolution ranged from 0.1-0.3 km s-1 for the JCMT data and were 0.05 and 0.1 km s-1 for the observations of respectively the C17O J=1-0 and J=2-1transitions at the IRAM 30 m. The line properties are summarized later in Sect. 5.
|Y||Ratio of the outer (r2) to inner (r1) radius|
|Optical depth at 100 m|
|Density power-law exponent|
|T1||Temperature at inner boundary (250 K)|
|Temperature of star (5000 K)|
Although these parameters may not seem the most straightforward choice, one of the advantages of DUSTY is the scale-free nature allowing the user to run a large sample of models and then compare a number of YSOs to these models just by scaling with distance and luminosity as discussed in Ivezic et al. (1999).
In determining the best fit model each of the calculated values are considered individually. The total
not make sense in a strictly statistical way, since the observations
going into these cases are not 100% independent. Another reason for
not combining the values of
into one total
is that the parameters
constrained by the SED and brightness profiles are different. For
example the brightness profiles provide good constraints on as seen in a 2D contour ()
plot of, e.g.,
Fig. 2, while these do not depend critically on the value
The most characteristic feature of the
-values for the SEDs on the other hand is the band of possible
models in contour plots for (
), giving an almost
one-to-one correspondence for a best fit
for each value
|Figure 2: contour plots for the modelling of L483-mm. In the four upper panels is fixed at 0.2, in the 4 middle panels is fixed at 0.9 and while in the 4 lower panels Y is fixed at 1400. The solid (dark) contours indicate the confidence limits corresponding to 1, 2etc.|
|Figure 3: The observed brightness profile for L483 at 450 m (upper panel) and 850 m (lower panel) with the best-fit models overplotted (full line). The dashed line indicates the beam profile used in the modelling.|
|Figure 4: The spectral energy distribution of L483: the points indicate the data, the line the best-fit model.|
These features are actually easily understood: and are the normalized profiles and should thus not depend directly on the value of . On the other hand since the peaks of the SEDs are typically found at wavelengths longer than 100 m, increasing and thus the flux at this wavelength, require the best-fit SED to shift towards 100 m, i.e., with less material in the outer cool parts of the envelope, which can be obtained by a steeper value of . The value of Y is less well constrained, mainly because of its relation to the temperature at the inner radius (and through that the luminosity of the central source). As illustrated in Fig. 2, Y is constrained within a factor 1.5-2.0 at the level, but the question is how physical the outer boundary of the envelope is: is a sharp outer boundary expected or rather a soft transition as the density and temperature in the envelope reaches that of the surrounding molecular cloud? In the first case, a clear drop of the observed brightness profile should be seen compared with a model with a (sufficiently) large value of Y, e.g., corresponding to an outer temperature of 5-10 K, less than the temperature of a typical molecular cloud. In the other case, however, such a model will be able to trace the brightness profile all the way down to the noise limit. The modelling of the CO lines (see Sect. 5) indicates that significant ambient cloud material is present towards most sources, so the transition from the isolated protostar to the parental cloud is likely to be more complex than described by a single power-law.
The features of the values of , and provide an obvious strategy for selecting the best fit models: first the best fit value of is selected on the basis of the brightness profiles and second the corresponding value of is selected from the contour plots. For a few sources there is not 100% overlap between the 450 and 850 m brightness profiles and it is not clear which brightness profile is better. The beam at 850 m is significantly larger than that at 450 m (15 vs. 9) and even though the beam is taken into account explicitly, the sensitivity of the 850 m data to variations in the density profile must be lower as is seen from Fig. 2. On the other hand, the 450 m data generally suffer from higher noise, so especially the weak emission from the envelope, which provides the better constraints on the outer parts of the envelope and thus the power-law exponent, will be more doubtful at this wavelength. The power-law slopes found from modelling the two brightness profiles agree, however, within the uncertainties ( ).
|Figure 5: Composite showing the brightness profiles of the sources overplotted with the best fit model at 450 m. The dashed line indicate the beam profile.|
|Figure 6: As in Fig. 5 but for the 850 m data.|
|Figure 7: Composite showing the SEDs of each source overplotted with the SED from the best fit model. The individual SEDs are based on literature searches, with the main references being Shirley et al. (2000) (class 0 objects and pre-stellar cores), Chandler & Richer (2000) (NGC1333-I2), Sandell et al. (1991) (NGC1333-I4A,B; 350 and 800 m), IRAS Faint Source Catalog (TMR1, L1489) and 1.3 mm fluxes from Motte & André (2001) for sources included in their sample.|
The derived power-law indices are for most sources in agreement with the predictions from the inside-out collapse model of . A few YSOs (esp. L1527, L483, CB244 and L1448-I2) have density distributions that are flatter, but as argued in Sect. 4.2, this can be explained by the fact that these sources seem to have significant departures from spherical symmetry. Apart from these sources, there is no apparent distinction between the class 0 and class I sources in the sample as shown in the plot of power-law slope versus bolometric temperature in the upper panel of Fig. 8. There is, however, a significant difference in the masses derived for these types of objects, with the class 0 objects having significantly more massive envelopes (lower panel of Fig. 8). This is to be expected since the bolometric temperature measures the redness (or coolness) of the SED, i.e., the amount of envelope material. Finally, there is no dependence on either mass or power-law slope with distance, which strengthens the validity of the derived parameters.
|Figure 8: The slope of the density distribution (upper panel) and the mass within the 10 K radius (lower panel) vs. the bolometric temperature of the sources. Class 0 objects are marked by " '', class I objects by " '' and pre-stellar cores by " ''. VLA1623 and IRAS 16293-2422 have been singled out with respectively "'' and "''. The mass of the pre-stellar cores in the lower panel is the mass of the Bonnor-Ebert sphere adopted for the line modelling.|
Chandler & Richer (2000) assumed that the envelopes were optically thin. In
this case, the temperature profile can be shown to be a simple
power-law as well and Chandler & Richer subsequently derived
analytical models which could be fitted directly to the SCUBA
data. For the sources included in both samples the derived power-law
indices agree within the uncertainties. However, as illustrated in
Fig. 9, the optically thin assumption for the
temperature distribution is not valid in the inner parts of the
envelope, especially of the more massive class 0 sources, so actual
radiative transfer modelling is needed to establish the temperature
profile, crucial for calculations of the molecular excitation and
|Figure 9: The temperature profiles for four selected sources with temperature profiles calculated in the Rayleigh-Jeans limit for the optically thin assumption overplotted for different dust opacity laws, . The dashed line corresponds to and the dash-dotted line to .|
Disk emission can contribute to the fluxes of the innermost points on the brightness profiles and thus lead to a steeper density profile. This is likely to be more important in the sources with the less massive envelope, i.e. class I objects. In tests where the fluxes within the innermost 15 of the brightness profiles are reduced by 50%, the best fit values of are reduced by 0.1-0.2. This is comparable to the uncertainties in the derived value of , but can introduce a systematic error.
It is interesting to note that there is no clear trend in the slope of the density profile with type of object. In the framework of the inside-out collapse model, one would expect a flattening of the density profiles, approaching 1.5 as the entire envelope undergoes the collapse. This is not seen in the data - actually the average density profile for the class I objects is slightly steeper than for the class 0 objects. On the other hand, the suggestion of an outer envelope with a flat density distribution and a significant fraction of material as suggested by Jayawardhana et al. (2001) can also not be confirmed by this modelling. As seen from the fits in Figs. 5 and 6 the brightness profiles only suggest departures from the single power-law fits in the outer regions in a few cases, so if such a component is present, it is not traced directly by the SCUBA maps. The slightly flatter density profiles for the class 0 objects could be a manifestation of such an outer component - but the density distributions of the sources in this sample are typically much steeper (i.e. ) than those modelled by Jayawardhana et al. ( ).
Myers et al. (1998) investigated the results of departure from spherical symmetry of an envelope when calculating the bolometric temperature of YSOs seen under various inclination angles. With a cavity in the pole region, roughly corresponding to the effect of a bipolar outflow, they found that the bolometric temperature could increase by a factor 1.3-2.5 for a typical opening angle of 25. A similar line of thought can be applied to our modelling: departure from spherical symmetry by having a thinner polar region will affect the determination of the SED - a source viewed more pole-on will see warmer material which leads to an SED shifted towards shorter wavelengths and accordingly a lower value of . Myers et al. also argued that the effect on a statistical sample would be rather small, e.g., compared to differences in optical depth, but for studies of individual sources like in our case, this effect might be of importance.
The brightness profile will also change in the aspherical case. In the case of a source viewed edge-on this would result in elliptically shaped SCUBA images with the 850 m data showing a more elongated structure since the smaller optical depths material in the polar regions would reveal material being warmer and thus having stronger emission at 450 m, compensating for the lack of material in these images. If we consider the case where the source is viewed entirely pole-on, the image would still appear circular, but the brightness profiles would show a steeper increase towards the center in the 450 m data (closer to the spherical case) than the 850 m data for the same reasons as mentioned above.
The good fits to the brightness profiles given the often non-circular
nature of the SCUBA images is expected based on the mathematical
nature of power-law profiles. Consider as an example a 2D image of a
source described by:
As a test, brightness profiles were extracted in angles covering respectively the flattest and steepest direction of L483 and L1527. Modelling the so-derived brightness profiles gives best-fit density profiles of 0.9 and 1.2, compared to the 0.9 derived as the average over the entire image for L483, and 0.6 and 0.8 for L1527 compared to 0.6 from the entire image. Thus, the brightness profile and the derived density profile might indeed be flattened by the asymmetry and could account for the somewhat flatter profiles found towards some sources. The discrepancies between the profiles along different directions are, however, not much larger than the uncertainties in the derived power-law slope, so these sources could have intrinsic flatter density distributions instead; with the present quality of the data both interpretations are possible.
Evans et al. calculated temperature distributions of Bonnor-Ebert spheres with radiative transport calculations, finding that the dust temperatures decline from about 13 K on the outside to about 7 K at the center. Ward-Thompson et al. (2002) examined ISOPHOT 200 m data for a sample of pre-stellar cores (including L1689B and L1544) and found that none of the cores had a central peak in temperature and that they could all be interpreted as being isothermal or having a temperature gradient with a cold center as result of external heating by the interstellar radiation field. Zucconi et al. (2001) derived analytical formulae for the dust temperature distributions in pre-stellar cores showing that these cores should have temperatures varying from 8 K in the center to around 15 K at the boundary. These equations will be useful for more detailed modelling of the continuum and line data, but for the present purpose the isothermal models are sufficient.
Modelling the pre-stellar cores using our method is not possible as is illustrated by the best fits for the pre-stellar cores shown in Figs. 5-7. Given the observational evidence that these cores do not have central source of heating, what is implicitly assumed in the DUSTY modelling, it is on the other hand comforting that our method indeed distinguishes between these pre-stellar cores and the class 0/I sources and that the obtained fits to the brightness profiles of the latter sources are not just the results of, e.g., the convolution of the "real'' brightness profiles with the SCUBA beam.
For modelling of sources without central heating Evans et al. note that their modelling does not rule out a power-law envelope density distribution for L1544, as opposed to, e.g., the case for L1689B: this would imply an evolutionary trend of the pre-stellar cores having a Bonnor-Ebert density distribution, which would then evolve towards a power-law density distribution with an increasing slope as the collapse progresses. In the modelling of spectral line data for the pre-stellar cores we adopt an isothermal Bonnor-Ebert sphere with for both L1689B and L1544 as this was the best fitted model in the work of Evans et al. (2001).
The 1D Monte Carlo code developed by Hogerheijde & van der Tak (2000) is
used together with the revised collisional rate coefficients from
Flower (2001), and a constant fractional abundance over the entire
range of the envelope is assumed as a first approximation. Furthermore
the dust and gas temperatures are assumed identical over the entire
envelope and any systematic velocity field is neglected. In the outer
parts of the envelopes, the coupling between gas and dust may break
down (e.g. Ceccarelli et al. 1996; Doty & Neufeld 1997) leading to differences
between gas and dust temperatures of up to a factor of 2.
|(km s-1)(a)||(km s-1)||Obs(b)||Mod(c)||Obs(b)||Mod(c)|
With the given physical model and assumed molecular properties, there
are two free parameters which can be adjusted by minimizing to model the line profiles for each molecular transition: the
fractional abundance [X/H2] and the turbulent line width .
Since emission from the ambient molecular cloud might contribute
to the lower lying 1-0 transitions, the derived abundances are based
only on fits to the 2-1 and 3-2 lines. The observed and modelled
line intensities are summarized in
Tables 7-10, whereas the parameters
for the best fit models are given in Table 11. While a
turbulent line width of 0.5-1.0 km s-1 is needed for most
sources to fit the actual line profile, the modelled line strengths
are only weakly dependent on this parameter compared to the fractional
abundance. For example, for L723 one derives C18O
km s-1 (
km s-1) to
km s-1), which
should be compared to the [C18O/H2]
km s-1 (
km s-1) quoted in
Table 11. The C18O and C17O spectra for
L723 and N1333-I2 with the best fit models
overplotted are shown in Fig. 10. The revised
collisional rate coefficients from Flower (2001) adopted for this
modelling typically change the derived abundances by less than 5%
compared to results obtained from simulations with previously
published molecular data. The relative intensities of the various
transitions and so the quality of the fits remain the same.
|Figure 10: The C18O and C17O spectra for L723 (upper four panels) and N1333-I2 (lower four panels) overplotted with the best fit models from the Monte Carlo modelling.|
As it is evident from Tables 7-10, the 2-1 and 3-2 lines can be well modelled using the above approach, while the 1-0 lines are significantly underestimated from the modelling, especially in the larger Onsala 20 m beam. This indicates that the molecular cloud material may contribute significantly to the observations or that the assumed outer radius is too small. The importance of the latter effect was tested by increasing the outer radius of up to a factor 2.5 for a few sources and it was found that while the 1-0 and 2-1 line intensities could increase in some cases by up to a factor of 2, the 3-2 line intensities varied by less than 10%, illustrating that the 3-2 line mainly trace the warmer (30 K) envelope material.
The importance of the size of the inner radius was tested as well. For L723 fixing the inner radius at 50 AU rather than 8 AU increases the best fit abundance by 5% to [C18O/H2] without changing the quality of the fit significantly (drop in of 0.1). This is also found for other sources in the sample and simply illustrates that increasing the inner radius corresponds to a (small) decrease in mass of the envelope, so that a higher CO abundance is required to give the same CO intensities. Yet, it is good to keep in mind that the inner radius of the envelope may be different, and even though its value does not change the results for CO, it might for other molecules, e.g., CH3OH and H2CO, which trace the inner warm and dense region of the envelope.
|Source||C17O 1-0||C17O 2-1|
|Source||C18O 1-0||C17O 1-0|
|Source||C18O 3-2(a)||C17O 3-2(a)|
The derived abundances are uncertain due to several factors. First, the physical model and its simplicity and flaws as discussed above. Second, the assumed dust opacities will affect the results: the dust opacities may change with varying density and temperature in the envelope, and will depend on the amount of coagulation and formation of ice mantles. These effects tend to increase with increasing densities or lower temperatures, which will lower the derived mass and thus increase the abundances necessary to reproduce the same line intensities. Comparing the models of dust opacities in environments of different densities and types of ice mantles as given by Ossenkopf & Henning (1994) indicates, however, that such variations should be less than a factor of two and so cannot explain the differences in the derived CO abundances.
The unknown contribution to the lower lines from e.g. the ambient molecular cloud may affect the interpretation. Even for the 2-1 transitions the cloud material may contribute, leading to higher line intensities than predicted from the modelling. Indeed as judged from Tables 7-10 there is a trend that while the intensities of 2-1 lines are underestimated by the modelling, the 3-2 lines are overestimated. This then means that the derived CO abundances are upper limits - at least for the warmer regions traced by the 3-2 lines.
As noted above the gas temperature could be lower by up to a factor of two in the outer parts of the envelope (Ceccarelli et al. 1996; Doty & Neufeld 1997), but this again would mainly affect the lines tracing the outer cold regions i.e. the 1-0 and 2-1 lines compared to the 3-2 lines. At the same time the gas temperatures from these detailed models may not be appropriate for our sample. Our sources all have lower luminosities than modelled e.g. by Ceccarelli et al. (1996), leading to envelopes that are colder on average, so that the relative effects of the decoupling of the gas and dust are smaller. Another effect is external heating, which as discussed in Sect. 4.3 may lead to an increase in the dust and gas temperatures in the outer parts.
A point of concern is also whether steeper density gradients could change the inferred abundances in light of the discussion about geometrical effects (Sect. 4.2). Steepening of the density distribution would tend to shift material closer to the center and so towards higher temperatures. This would correspondingly change the ratio between the 2-1 and 3-2 lines towards lower values (stronger 3-2 lines) and more so for the C18O data than the C17O data, since it traces the outer less dense parts of the envelope. Altogether none of the effects considered change the conclusion that the abundances in the class 0 objects are lower than those found for the class I objects.
One important conclusion of the derived abundances is that one should
be careful when using CO isotopes to derive the H2 mass of, e.g.,
envelopes around young stars assuming a standard abundance. With
depletion the derived H2 envelope masses will be
underestimated. Another often encountered assumption, which may
introduce systematic errors, is that the lower levels are thermalised
and that the Boltzmann distribution can be used to calculate the
excitation and thus column density of a given molecular species. In
Fig. 12, the level populations for C18O in the outer shell
of the model of two sources, L723 and TMR1, are
shown together with the variation of the ratios of the level
populations from the Monte Carlo modelling by assuming LTE. For the
more dense envelope around the typical class 0 object, L723,
the LTE assumption gives accurate results within 5-10% in the lower
levels, but somewhat more uncertain results for the higher levels. For
the less dense envelope around TMR1 (class I object) it is,
however, clear that the levels are subthermally excited and the LTE
approximation provides a poor representation of the envelope
|Figure 11: The fitted abundances vs. envelope mass (left) and bolometric temperature (right) of each source for respectively the C18O data (top) and C17O data (bottom). The sources have been split into groups (class 0, class I and pre-stellar cores) with VLA 1623 and IRAS 16293-2422 separated out using the same symbols as in Fig. 8: class 0 objects are marked with " '', class I objects with " '', pre-stellar cores with " '', VLA 1623 with "'' and IRAS 16293-2422 with "''. The vertical lines in the figures illustrate the abundances in quiescent dark clouds from Frerking et al. (1982).|
|Figure 12: The derived level populations for C18O in the envelopes around the class 0 object L723 (left column) and class I object TMR1 (right column). In the upper panels the level populations from the modelling in the outer shell of the envelope are shown (dashed line) together with the predictions from the LTE assumptions (solid line). In the lower panels, the ratio of the resulting level populations from the modelling and the LTE predictions are shown with varying density for the J=0 to 3 levels.|
The apparently low CO abundances and the possible relation to freeze-out of CO raises the question whether the assumption of a constant fractional abundance is realistic: freezing out of pure CO-ice and isotopes is expected to occur at roughly 20 K under interstellar conditions (e.g., Sandford & Allamandola 1993), so one would expect to find a drastic drop in CO abundance in the outer parts of the envelope. Due to the uncertainties in the properties of the exterior regions, however, a change in abundance at 20 K can neither be confirmed nor ruled out. As Table 7 indicates both the intensities of the 2-1 and 3-2 lines of C18O and C17O can be fitted with a constant fractional abundance for most sources. One can introduce strong depletion through a "jump'' in the fractional abundance in the region of the envelope with temperatures lower than 20 K. This naturally leads to lower line intensities of especially the 2-1 and 1-0 lines, but also the 3-2 lines. Since it is mainly the 3-2 line that constrains the abundance in the inner part, it is possible to raise the abundances of the warm regions by up to a factor 2, if CO in the outer part is depleted by a factor of 10 or more. The modelled 1-0 and 2-1 line intensities then, however, also become weaker, which one has to compensate for by introducing even more cold (depleted) material outside the 10 K boundary, accounting for 50% or more of the observed 2-1 line emission and almost all the 1-0 emission.
If the dust opacity law varies with radial distance in the envelope, increasing with the higher densities, the CO abundances in the warmer regions could be higher by a factor of two. Combining these effects would then lower the derived CO evaporation temperature in the envelope material towards the 20 K evaporation temperature of pure CO ice. A changing dust opacity law would be an interesting result in itself, which possibly can be confirmed or disproved through modelling of other chemical species and/or by comparing the exact line profiles with realistic models of the velocity field. Without such an effect introduced, however, the results presented here favors the somewhat higher evaporation temperature for the CO.
The differences between the class 0 and I objects and a warmer evaporation temperature may be consistent with new laboratory experiments on the trapping and evaporation of CO on H2O ice by Collings et al. (2002). Their experiments refer to amorphous H2O ice accreted layer-by-layer so that it has a porous ice structure. This situation may be representative of the growth of H2O ice layers in pre-stellar cores and YSO envelopes. CO is deposited on top of the H2O ice. When the sample is heated from 10 to 30 K (laboratory temperatures), some of the CO evaporates, but another fraction diffuses into the H2O ice pores. Heating to 30-70 K allows some CO to desorb, but a restructuring of the H2O ice seals off the pores, and the remaining CO stays trapped until at least 140 K. Under interstellar conditions, the temperatures for these processes may be somewhat lower, but it does indicate that a significant fraction of CO can be trapped in a porous surface and that evaporation may occur more gradually from 30 K to 60 K. A similar property was suggested by Ceccarelli et al. (2001), who found that the dust mantles in the envelope around IRAS 16293-2422 had an onion-like structure with H2CO being trapped in CO rich ices in the outermost regions and with the ices becoming increasingly more H2O rich when moving inwards toward higher temperatures. In this scenario, it is not surprising that even the 3-2 lines tracing the warmer material indicate low CO abundances. Observations of even higher lying CO rotational lines (e.g., 6-5) would be needed probe the full extent of this evaporation. Also infrared spectroscopy of CO ices may reveal differences between the class 0 and I objects.
The physical models derived using this method have been applied in Monte Carlo modelling of C18O and C17O data, adopting an isothermal Bonnor-Ebert sphere as a physical model for the pre-stellar cores. The 2-1 and 3-2 lines can be modelled for all sources with constant fractional abundances of the isotopes with respect to H2and an isotope ratio [18O/17O] of 3.9, in agreement with the "standard'' value for the local interstellar medium. The 1-0 lines intensities, however, are significantly underestimated in the models compared to the observations, indicating that ambient cloud emission contributes significantly or that the outer parts of the envelopes are not well accounted for by the models. The derived abundances increase with decreasing envelope mass - with an average CO abundance of for the class 0 objects and pre-stellar cores, and for the class I objects. The 3-2 lines indicate that the lower CO abundance in class 0 objects also applies to the regions of the envelopes with temperatures higher than 20-25 K, the freeze-out of pure CO ice. This feature can be explained if a significant fraction of the solid CO is bound in a (porous) ice mixture from where it does not readily evaporate. The physical models presented here will form the basis for further chemical modelling of these sources.
The authors thank Michiel Hogerheijde and Floris van der Tak for use of their Monte Carlo code and stimulating discussions, and Helen Fraser for sharing the results of the new experiments on CO evaporation prior to publication. They are grateful to Sebastien Maret for carrying out the IRAM observations in November 2001. The referee, Cecilia Ceccarelli, is thanked for thorough and insightful comments, which helped in clarifying the results of the paper. The work of JKJ is funded by the Netherlands Research School for Astronomy (NOVA) through a netwerk 2, Ph.D. stipend. FLS is funded by grants from the University of Leiden and the Netherlands Foundation for Scientific Research (NWO) no. 614.041.004. Astrochemistry in Leiden is supported by a Spinoza grant. This article made use of data obtained through the JCMT archive as Guest User at the Canadian Astronomy Data Center, which is operated by the Dominion Astrophysical Observatory for the National Research Council of Canada's Herzberg Institute of Astrophysics. The SIMBAD database, operated at CDS, Strasbourg, France has been intensively used as well.