A&A 518, A52 (2010)

Chemical study of intermediate-mass (IM) Class 0 protostars

CO depletion and N₂H⁺ deuteration

T. Alonso-Albi¹ - A. Fuente¹ - N. Crimier² - P. Caselli³ - C. Ceccarelli² - D. Johnstone^4,5 - P. Planesas^1,6 - J. R. Rizzo⁷ - F. Wyrowski⁸ - M. Tafalla¹ - B. Lefloch² - S. Maret² - C. Dominik⁹

1 - Observatorio Astronómico Nacional (OAN, IGN), Apdo 112, 28803 Alcalá de Henares, Spain
2 - Laboratoire d'Astrophysique, Observatoire de Grenoble, BP 53, 38041 Grenoble Cedex 9, France
3 - School of Physics & Astronomy, E.C. Stoner Building, The University of Leeds, Leeds LS2 9JT, UK
4 - Department of Physics & Astronomy, University of Victoria, Victoria, BC, V8P 1A1, Canada
5 - National Research Council of Canada, Herzberg Institute of Astrophysics, 5071 West Saanich Road, Victoria, BC, V9E 2E7, Canada
6 - Joint ALMA Observatory, El Golf 40, Las Condes, Santiago, Chile
7 - Centro de Astrobiología (CSIC/INTA), Laboratory of Molecular Astrophysics, Ctra. Ajalvir km. 4, 28850, Torrejón de Ardoz, Spain
8 - Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
9 - Anton Pannekoek Astronomical Institute, University of Amsterdam, PO Box 94249, 1090 GE Amsterdam, The Netherlands

Received 24 February 2010 / Accepted 16 April 2010

Abstract
Aims. We are carrying out a physical and chemical study of the protostellar envelopes in a representative sample of IM Class 0 protostars. In our first paper we determined the physical structure (density-temperature radial profiles) of the protostellar envelopes. Here, we study the CO depletion and N₂H⁺ deuteration.
Methods. We observed the millimeter lines of C¹⁸O, C¹⁷O, N₂H⁺ and N₂D⁺ towards the protostars using the IRAM 30m telescope. Based on these observations, we derived the C¹⁸O, N₂H⁺ and N₂D⁺ radial abundance profiles across their envelopes using a radiative transfer code. In addition, we modeled the chemistry of the protostellar envelopes.
Results. All the C¹⁸O 1 0 maps are well fit when assuming that the C¹⁸O abundance decreases inwards within the protostellar envelope until the gas and dust reach the CO evaporation temperature, 20-25 K, where the CO is released back to the gas phase. The N₂H⁺ deuterium fractionation in Class 0 IMs is [N₂D⁺]/[N₂H⁺]=0.005-0.014, two orders of magnitude higher than the elemental [D/H] value in the interstellar medium, but a factor of 10 lower than in prestellar clumps. Chemical models account for the C¹⁸O and N₂H⁺ observations if we assume the CO abundance is a factor of 2 lower than the canonical value in the inner envelope. This could be the consequence of the CO being converted into CH₃OH on the grain surfaces prior to the evaporation and/or the photodissociation of CO by the stellar UV radiation. The deuterium fractionation is not fitted by chemical models. This discrepancy is very likely caused by the simplicity of our model that assumes spherical geometry and neglects important phenomena like the effect of bipolar outflows and UV radiation from the star. More important, the deuterium fractionation is dependent on the ortho-to-para H₂ ratio, which is not likely to reach the steady-state value in the dynamical time scales of these protostars.

Key words: ISM: abundances - ISM: clouds - stars: formation - circumstellar matter

1 Introduction

Intermediate-mass young stellar objects (IMs) share many characteristics with high-mass stars (clustering, PDRs) but their study presents an important advantage: many are located closer to the Sun ( kpc) and in less complex regions than massive star-forming regions. On the other hand, they are also important for understanding planet formation since Herbig Ae stars are the precursors of Vega-type systems. Despite this, IMs have been studied very little so far. A few works on HAEBE stars have been carried out at millimeter wavelengths (Fuente et al. 1998, 2002; Henning et al. 1998), but almost nothing has been done for their precursors, the Class 0 IM objects.

Chemistry has been successfully used to determine the physical structure and investigate the formation and evolution of low-mass YSOs. Chemical diagnostics have also been shown to be good indicators of the protostellar evolution in these objects (see e.g. Maret et al. 2004; Jørgensen et al. 2005). However, very few works deal with IMs. Fuente et al. (2005a) present a chemical study of the envelopes of the Class 0 IM protostar NGC 7129-FIRS 2 and the young Herbig Be star LkH 234. They find that the changes in the physical conditions of the envelope during its evolution from the Class 0 to the Class I stage (the envelope is dispersed and warmed up) strongly influence the molecular chemistry. The Class 0 object NGC 7129-FIRS 2 presented evidence of H¹³CO⁺ depletion. Moreover, the deuterium fractionation, measured as the DCO⁺/H¹³CO⁺ratio, decreases by a factor of 4 from the Class 0 to the Herbig Be star, very likely owing to the increase in the kinetic temperature. Regarding the abundance of complex molecules, the beam-averaged abundances of CH₃OH and H₂CO increase from the Class 0 to the Herbig Be star. A hot core was also detected in NGC 7129-FIRS 2 (Fuente et al. 2005a,b). Although two objects are not enough to establish firm conclusions, these pioneering results suggest that chemistry is also a good indicator of the evolution of IMs.

We are carrying out a chemical study of a representative sample of IM Class 0 YSOs This is the first systematic chemical study of IM Class 0 objects that has been carried out so far. Some properties like the temperature of the protostellar envelope and the clustering degree depend on the final stellar mass, so the results for low-mass stars cannot be directly extrapolated to intermediate-mass objects. In the first paper (Crimier et al. 2010, hereafter C10), we determined the physical structure (density-temperature radial profiles) by modeling the dust continuum emission. We now present the observations of the millimeter lines of C¹⁸O, C¹⁷O, N₂H⁺, and N₂D⁺ in the same sample. Our goal is to investigate the CO depletion and N₂H⁺ deuteration in these Class 0 YSOs. For comparison, we also include 2 Class I objects, LkH 234, and S140.

Table 1:   Selected sample.

2 Observational strategy

Our selection was made to have a representative sample of Class 0 IM YSOs, including targets with different luminosities (40-10 ) and evolutionary stages. An important complication in the study of massive stars is that they are located in complex regions and are therefore difficult to model. The targets in this sample were chosen to lie preferentially in isolated areas with respect to the 30 m telescope beam. We also selected sources for which continuum maps at submillimeter and/or millimeter wavelengths are available in order to be able to model their envelopes. The list of sources and their coordinates are shown in Table 1. To provide a comparison with Class I sources, we added S140 and LkH 234 to the sample.

Most of the observations reported here were carried out with the IRAM 30 m telescope at Pico de Veleta (Spain) during three different observing periods in June 2004 in position switching mode. Our strategy was to first make long integration single-pointing observations towards the star position and then to make 96 maps around the center position. The maps were sampled with a spacing of 12'' in the inner 48'' regions and 24'' outside. The only exceptions were L1641 S3 MMS1 and S140. In L1641 S3 MMS1, we only observed a radial strip in C¹⁷O and C¹⁸O. We did not observe C¹⁷O and C¹⁸O maps in S140. A summary of the observations for each source is shown in Table 1, and the list of observed lines and the telescope characteristics is shown in Table 2. During the observations, lines of the same species were observed simultaneously using the multireceiver capability of the 30 m telescope. In this way, we minimized relative pointing and calibration errors. As backends we used in parallel an autocorrelator split into several parts providing a spectral resolution that was always better than 78 kHz and a 1 MHz-channel-filter-bank. Examples of the single-pointing observations are shown in Figs. 1 and 2. The intensity scale is the main brightness temperature.

Observations of the N₂H⁺ line (  GHz) were carried out towards Cep E-mm, IC 1396 N, NGC 7129-FIRS 2, and L1641 S3 MMS1 using the JCMT telescope at the Mauna Kea (Hawaii). In all these sources, we carried out small maps of 75 with a spacing of 25'', but the emission was only detected towards the star position (see Fig. 3). All the lines were observed with a spectral resolution of 0.488 MHz. The spectra of the line are shown in Fig. 3.

In this paper, we also use the C¹⁸O, N₂H⁺ and N₂D⁺ data towards NGC 7129-FIRS 2 and LkH 234, which has already been published by Fuente et al. (2005b).

Table 2:   Description of the observations.

3 Results

The spectra of the C¹⁸O 1 0, C¹⁷O 1 0, N₂H⁺ 1 0 and N₂D⁺ 2 1 lines towards the star position are shown in Figs. 1 and 2. The integrated intensity maps of the C¹⁸O 1 0, N₂H⁺ 1 0, and N₂D⁺ 2 1 lines are shown in Figs. 4 (Class 0) and 5 (Class I). In all the Class 0 sources, the emission of the N₂H⁺ 1 0 line is compact and peaks towards the star position, revealing a dense core around the star. However, towards the Class I sources, the peaks of the N₂H⁺ emission are offset from the far-IR source. For the Herbig Be star LkH 234, the emission peak of the N₂H⁺ 1 0 line is located at an offset (-6'',18'') from the star position (see Fuente et al. 2005a and Fig. 5). In S140, the N₂H⁺ emission peaks in an arc-shaped feature surrounding the sources IRS 1, 2 and 3. In Class I sources, either the N₂H⁺ is destroyed by the evaporated CO or the outflow, or the UV radiation from the star has already disrupted the parent core.

The emission in the C¹⁸O 1 0 line usually presents a different morphology from that of N₂H⁺ 1 0 line. In CB3 and IC 1396N, the emission in the C¹⁸O  line has an elongated shape, much more extended than that of N₂H⁺ (see Fig. 4). In OMC2 FIR 4 and NGC 7129-FIRS 2, the emission of the C¹⁸O line surrounds the star position instead of having a maximum towards it (see also Fuente et al. 2005a). Only in the low-luminosity sources Serpens-FIRS 1 and Cep E-mm does the emission from the C¹⁷O and C¹⁸O lines peak at the star position and present a morphology similar to that of the N₂H⁺ emission.

For the sources in which the signal-to-noise ratio of the N₂D⁺ map is high enough, Serpens-FIRS 1 and NGC 7129 FIRS 2, we compared the N₂D⁺ 2 1 and N₂H⁺ 1 0 maps. In Serpens-FIRS 1, the N₂D⁺ 2 1 emission does not peak towards the star, but does in the case of NGC 7129-FIRS 2.

In Fig. 6, we show the mean integrated intensity emission of the C¹⁸O 1 0, N₂H⁺ 1 0, and N₂D⁺ 2 1 lines in concentric rings around the protostar. In all the sources, the emission of the N₂H⁺ 1 0 line decreases outwards from the star. The radial profiles of the C¹⁸O 1 0 line, however, greatly differ from one source to the next. The presumably youngest sources, OMC2 FIR 4 and NGC 7129-FIRS 2, show a flat profile. This is the expected picture when the CO abundance decreases towards the center, mainly because the molecules are frozen onto the grain surfaces. In Serpens-FIRS 1, Cep E-mm, IC 1396 N, and CB3, the C¹⁸O 1 0 emission decrease with the distance from the star, similar to N₂H⁺ but with a less steep profile. Intense emission of the C¹⁸O 1 0 line is detected at the edges of the protostar envelopes, which shows that a lower density envelope also contributes to the emission of this line. This envelope contribution is especially important for Serpens-FIRS 1. In fact, we must model the envelope contribution in order to fit the C¹⁸O 1 0 emission from this protostellar core (see Appendix A1). One could interpret the different radial profiles of the C¹⁸O 1 0 emission as an evolutionary trend, with the youngest protostars having flat profiles while the oldest have steep profiles. Given the complexity of these regions, however, one should be cautious and use multiple evolutionary tracers to define relative age. For instance, a lack of C¹⁸O emission towards the star position could come from CO depletion in the case of a very young object or to photodissociation for a borderline Class 0/I. The radial profiles of the N₂D⁺ 2 1 (3 2 for NGC 7129-FIRS 2) are more uncertain because of the low S/N ratio of the maps.

4 LTE column densities

We have derived the C¹⁸O column densities using the rotation diagram method. This method gives an accurate estimate of the column density provided that the emission is optically thin, is thermalized, and arises from a homogeneous and isothermal slab. In the case of a density and temperature distribution, the derived column density represents an average value over the observational beam.

The excitation of C¹⁸O and C¹⁷O are very similar, and we have better signal-to-noise spectra for the more intense C¹⁸O lines. For this reason, we derived the C¹⁷O column density by assuming optically thin emission and the rotation temperature derived from the C¹⁸O data. The molecule N₂H + presents hyperfine splitting. This allows us to estimate the line opacity directly from the hyperfine line ratios. We derived the total N₂H⁺ column density from the opacity of the N₂H⁺ 1 0 line, and assumed the rotation temperatures derived from N₂D⁺data when there was no N₂H⁺ 4 3 observations. The same rotation temperature was always used for both N₂H⁺ and N₂D⁺. We assumed a beam filling factor of 1 for all the lines regardless of the observational beam size. This approach overestimates the rotation temperature when the emission is centrally peaked. In Table 3 we present the derived C¹⁸O, C¹⁷O, N₂H⁺, and N₂D⁺ column densities. We also show the N(C¹⁸O)/N(C¹⁷O) (hereafter R₁) and the N(N₂D⁺)/N(N₂H⁺) (hereafter R₂) ratios. For the whole sample, we obtain a value for R₁ around 3.3+/-0.3, the expected value for optically thin emission. In the worst case, R₁=2.5, the correction to the C¹⁸O column density because of the opacity is only a factor 1.3. Thus, opacity effects are not important in our C¹⁸O column density estimates.

Table 3:   LTE column densities.

4.1 N₂D⁺ and N₂H⁺

The deuterium fractionation of N₂H⁺ (R₂) strongly depends on the CO depletion factor and the gas temperature (Caselli et al. 2002; Ceccarelli & Dominik 2005; Daniel et al. 2007). In our sample of IM Class 0 protostars, we derive values of R₂ ranging from 0.005 to 0.014. These values are 3 orders of magnitude higher than 10^-5, the elemental value in the interstellar medium (Oliveira et al. 2003). Nevertheless, the R₂ values are a factor of 10 lower than those found in prestellar clumps by Crapsi et al. (2005). According to the values of R₂, we can classify our sources into two groups: (i) highly deuterated sources that have values of R₂ >0.01 and (ii) moderately deuterated sources with R₂ <0.01. Cep E-mm, CB3, and NGC 7129-FIRS 2 belong to the first group. Assuming that the N(N₂D⁺)/N(N₂H⁺) ratio is a good gas temperature indicator, these protostars should be the coldest and very likely the youngest of our sample. We have to be cautious with CB3, however. Since this source is the most distant (d=2500 pc), our single-pointing observations trace a larger fraction of the envelope. The sources Serpens-FIRS 1, IC 1396 N, OMC2-FIR 4, and S140 belong to the second group. S140 is a Class I YSO and IC 1396 N is considered a borderline Class 0/I YSO. The results are thus consistent with our interpretation of the objects in this group as being more evolved. In Serpens-FIRS 1, however, the N₂D⁺ 1 0 emission comes mainly from the lower density envelope and the ratio cannot be considered a tracer of the evolutionary stage of the protostar. OMC2 FIR 4 is also difficult to classify. The deuterium fractionation suggests an evolved object but the morphology of the C¹⁸O map is more consistent with a young Class 0 star. The peculiarity of this protostellar envelope has already been pointed out by Crimier et al. (2009). They find that the envelope of OMC2 FIR 4 is peculiarly flat and warm with a radial density power-law index of 0.6.

We present the N₂H⁺ deuterium fractionation as a function of the N(C¹⁸O)/N(N₂H⁺) ratio in Fig. 7. An excellent correlation between these two quantities is found in prestellar cores with the deuterium fractionation decreasing with increasing N(C¹⁸O)/N(N₂H⁺) ratio (see Crapsi et al. 2005). This correlation is not valid for IM Class  0 protostars. As discussed in the following sections, this is due to the complexity of these intermediate mass star-forming regions. In prestellar cores the chemistry of these species is only driven by the CO depletion. In IMs, other phenomena like photodissociation and shocks are also playing important roles.

Table 4:   Temperature-density (T-n) profiles from Crimier et al. (2010).

Table 5:   Abundance profiles derived from the observations.

5 Radiative transfer model

We utilized a general ray-tracing radiative transfer code to derive the fractional abundance profiles of C¹⁸O, N₂H⁺ and N₂D⁺ across the envelopes. Assuming appropriate radial profiles for the temperature, density, molecular abundance, and turbulence velocity, this model calculated the brightness temperature distribution on the sky. The model map was then convolved with the telescope beam profile. The underlying source geometry was assumed to be a sphere, with the inner and outer radii, and temperature-density (T-n) profiles derived by C10. The size of the grid was set to  cells. The cells have different sizes along the line of sight to account for the different slopes in the density and temperature profiles. We used very small cells (100 AU) near the center, where the temperature and density gradients are highest. In the outer regions, the cells reach several thousand AU in size. The turbulent velocity was assumed to be fixed at 1.5 km s^-1, consistent with the linewidths of the observed lines. Finally, in each cell, the excitation temperature was calculated with the RADEX code (van der Tak et al. 2007), which uses the slab LVG approximation at each shell. We used the collisional rates provided by the LAMDA database (Schöier et al. 2005). For N₂D⁺, we used the same collisional rates as for N₂H⁺.

With these assumptions, we searched for the best fit on every source and molecule observed. The fitting process consisted of averaging the observed and modeled fluxes in concentric rings around the center position using GILDAS software. The center position was selected to be the continuum emission peak at 850 m. The radius of the first ring was set to HPBW/4, and incremented by the same amount in consecutive rings. As a first step we searched for the best fit using a constant abundance. Since this approach seldom produced good fits, we next searched for the best fit using radial power functions and step functions. These functions were not selected arbitrarily but were selected to mimic the predictions of chemical models. The step function accounts for the abrupt sublimation of the CO ices thanks to the increase in the dust temperature going inwards, whereas the power-law profile accounts for a smoother change in the abundance of the species. The angular resolution of our observations (16''-27'' depending on the frequency, i.e, 16 000-27 000 AU at a distance of 1000 pc) prevents us from tracing the inner region (a few 1000 AU) of the protostellar envelope. For this reason, we assume a constant abundance in the inner part.

We fit the integrated intensity maps of the C¹⁸O 1 0, N₂H⁺ 1 0, N₂H⁺ 4 3, and N₂D⁺ 2 1 lines, when possible. The maps of the C¹⁸O 2 1 and C¹⁷O 1 0 and 2 1 lines are of less quality so we did not use them in our fits. Some sources show deviations from the spherical symmetry because of the contribution of the surrounding molecular cloud (Serpens-FIRS 1, CB3) or because of the cavity excavated by the outflow (IC 1396 N). In these cases we masked part of the map to avoid this contamination in the rings. A short description of earch source and the individual details about the modeling are given in Appendix A. The areas masked are shown in Fig. 4, and the abundance profiles obtained from the model are listed in Table 5. We did not model the Class I sources S140 and LkH 234. The modeling of OMC2 FIR 4 will be the subject of a forthcoming paper.

6 Limitations of our model: CO evaporation temperature

The model used here is not a reliable predictor of the chemical and physical properties within the inner envelope for several reasons. First of all we only considered the low-J rotational molecular lines. The emission of these transitions arises mainly from the outer envelope. This fact, together with the limited spatial resolution of our single dish observations, prevent us from determining the variation in the abundance in the inner envelope. Additionally, in our modeling we use the n-T profiles derived by C10. These profiles are a reasonable approximation for the outer envelope but are relatively unconstrained in the inner region. Some parameters, like the inner radius of the envelope and the dust temperature at this radius, are thus poorly determined.

Another important source of uncertainty is the degeneracy of the solutions. For instance, in the case of C¹⁸O we can obtain similar results by varying the molecular abundance in the inner region (X₀) or the CO evaporation radius (R₀), as long as the line is optically thin and the number of molecules in the beam remains constant. Here we discuss the impact of this degeneracy in the derived CO evaporation radius.

We ran a grid of models for Cep E-mm, IC 1396 N, and CB3. These are the sources in which the ratio between the angular diameter of the envelope and the HPBW of the telescope is greatest, so we are more sensitive to the spatial variations of the chemical abundances. In all these models, we assume that the C¹⁸O abundance profile is X=X₀ for radii less than the evaporation radius, R₀, and ) for radii larger than R₀. For each value of X₀, we fit R₀ and . X₀ is the un-depleted C¹⁸O abundance, R₀ the evaporation radius, and measures the gradient in the C¹⁸O abundance due to depletion. The canonical abundance of C¹⁸O in principle depends on the Galactocentric distance (DGC). Following Wilson & Matteucci (1992), the CO abundance is given by

Following Wilson & Rood (1994), the oxygen isotope ratio ¹⁶O/¹⁸O depends on DGC according to the relationship ¹⁶O/¹⁸O = 58.8  DGC (kpc) + 37.1. At the distance of the Sun, 8.5 kpc, . This is an average value, so there may be local effects. In our models we varied X₀ from to , the range of X₀ values that we consider reasonable. For each value of X₀, we varied R₀ and to find the best fit. In Fig. 8 we show the results of our models for Cep E. The best fit is , R₀=2510 AU, with  K km s^-1. We consider that the models with rms values departing by a factor of less than 1.3 from the minimum value are still acceptable. Following this criterion, we have two different families of solutions. Assuming , the best fit is obtained with and -3100 AU. Assuming , we have reasonably good solutions with and -3200 AU. We have not found any acceptable solution for higher values of X₀. The evaporation radius of C¹⁸O is thus  AU. These radii correspond to dust temperatures of 20-25 K. Interestingly, we would have a better fit to the observations (  K km s^-1) if we allowed the C¹⁸O abundance in the inner region to fall below . This is the solution we show in Table 5. We do not reject this solution because there are several physical reasons that the CO abundance could be lower in the inner envelope. First of all, C¹⁸O could have been photodissociated in the regions very close to the recently formed star. Alternatively, part of the CO could have been transformed into CH₃OH on the icy mantles, and the CO abundance once evaporated would then be different from the initial value. However, given the large uncertainty in the physical structure of the inner parts of the envelope, we cannot draw any conclusion.

In Fig. 9 we show the results for IC 1396 N. The best fit is , R₀=2000 AU, and with rms = 0.36 K km s^-1. We consider that the models with values of rms < 0.5 are acceptable. Following this criterion, we only find solutions for and R₀ <3000 AU; i.e., the CO evaporation temperature should be >23 K. Again, the fit can be significantly improved if values of X₀ lower than are considered. This is similar to the case of Cep E. Finally, in Fig. 10, we show the same grid of models for CB3. The best fit is for , R₀=4000 AU, and with rms = 2.0 K km s^-1, which corresponds to a CO evaporation temperature of 30 K. However, this is not a good fit. In fact, we obtain a much better solution when assuming a step function and  K (see Table 5).

For Serpens-FIRS 1 and NGC 7129 FIRS2, the protostars are barely resolved by our single-dish observations. In addition, the contribution of the surrounding cloud is very significant. For these reasons, we did not carry out a study similar to what is described above. We looked only for abundance profiles that fit our observations and are consistent with the behavior predicted by chemical models (see Table 5). In these two cases, we can fit the observations when assuming that the CO is evaporated at temperatures around 20-25 K. We conclude, therefore, that our C¹⁸O observations towards the Class 0 stars are better fit assuming CO evaporation temperatures of 20-25 K that are consistent with the CO being bound in a CO-CO matrix.

The fits to the N₂H⁺ and N₂D⁺ emission profiles have been done by hand. For N₂H⁺ we tried two options:(i) constant N₂H⁺ abundance and (ii) an abundance profile with a radial variation similar to that of C¹⁸O. Option (ii) follows the expectation that the N₂H⁺ abundance is strongly dependent on the CO abundance. In the case of CB3 and NGC 7129 FIRS 2, we are able to fit the N₂H⁺ observations by assuming a constant abundance. For the rest of sources, the N₂H⁺ observations are better fit by assuming that the N₂H⁺ is depleted in the cold regions similarly to CO.

One expects an annular distribution for the N₂D⁺ abundance, with the maximum N₂D⁺ abundance in the region with the greatest CO depletion. The morphology observed in the integrated intensity map of the N₂D⁺ 2 1 line towards Serpens-FIRS 1 suggests that the N₂D⁺ emission is dominated by the foreground molecular cloud. For this reason we selected a one-step function with a constant abundance inside and outside the protostellar core. In the case of NGC 7129-FIRS 2, we detected a clump of N₂D⁺ 3 2 emission towards the star. However, interferometric observations show that on scales of a few 1000 AU, there is a lack of emission towards the star (Fuente et al. 2005a). We thus fit the N₂D⁺ emission assuming an annular abundance distribution.

7 Chemical model

The chemical composition was modeled with the simple chemical code originally described in Caselli et al. (2002), and updated following Caselli et al. (2008) with new measurements of the CO and N₂ binding energies (Collings et al. 2003; Öberg et al. 2005) and sticking coefficients (Bisschop et al. 2006), thermal desorption, and the detailed physical structure. The clouds are assumed to be spherically symmetric, with the density and temperature profiles derived from the dust continuum emission. An interpolation procedure has been included in the code in order to have smooth profiles. The chemical network contains the neutral species CO and N₂, which can freeze out onto dust grains and desorb owing to cosmic ray impulsive heating (as in Hasegawa & Herbst 1993) and by thermal evaporation (following Hasegawa et al. 1992). The initial abundances of CO and N₂ have been fixed to (Frerking et al. 1982) and , respectively. The abundances of molecular and atomic nitrogen are difficult to determine in dense cores, but recent works suggest low values in low-mass, star-forming regions: X(N) = n(N)/n(H (Hily-Blant et al. 2010), and X(N (Maret et al. 2006). Our adopted value is consistent with 25% of the total abundance of N (n(N)/n(H , see Anders & Grevesse 1989) locked in N₂ and with the pseudo-time dependent models of Lee et al. (1996), after the chemistry reaches steady state. Although atomic oxygen can affect the amount of deuterium fractionation (see discussion in Caselli et al. 2002), no atomic oxygen is included in the code because of the large uncertainties associated with its value (see e.g. Caux et al. 1999 and Melnick & Bergin 2005). This issue is discussed further at the end of this section.

The abundances of the molecular ions (HCO⁺, N₂H⁺, H₃⁺, and all their deuterated forms) were calculated in terms of the instantaneous abundances of neutral species, assuming that the timescale for ion chemistry is much shorter than for freeze-out (Caselli et al. 2002). The rate coefficients are adopted from the UMIST database (http://www.udfa.net). For the proton-deuteron exchange reactions (such as H₃⁺ + HD H₂D⁺ + H₂), we used the rates measured by Gerlich et al. (2002), which better fit the deuterium fractionation in low-mass Class 0 sources, as recently found by Emprechtinger et al. (2009). Hugo et al. (2009) have recently measured the proton-deuteron rate coefficients again, finding average values of total rates (for H₃⁺ and its deuterated isotopologues, the total rate refers to the average of multiple rates weighted according to the fraction of ortho and para forms; Sipilä et al. 2010) about 4-5 times more than those derived by Gerlich et al. (2002; see also Sipilä et al. 2010). Although this difference could affect our results by increasing the deuterium fractionation by a similar factor (thus worsening the comparison with observations), we decided not to include the new values because the nuclear spin variants of all H₃⁺ isotopologues and H₂ have not been distinguished in the current model. In fact, as Pagani et al. (1992) and Flower et al. (2004) showed, small temperature variations significantly alter the ortho-to-para ratio of H₂, which in turn strongly affects the D-fractionation (ortho-H₂ being more efficient than para-H₂ in driving back the proton-deuteron exchange reactions thus decreasing the D-fractionation), even in the temperature regime between 9 and 20 K, where CO is mostly frozen onto dust grains (see also the discussion about the drop in ortho-H₂D⁺ column density at temperatures above 10 K in Caselli et al. 2008). Such temperature variations are definitely present in the envelopes of intermediate-mass Class 0 protostars (see also Emprechtinger et al. 2009), so that a simple increase in the rate coefficients without accounting for nuclear spin variants and, in particular, the possible increase in the H₂ ortho-to-para ratio may overestimate the D-fractionation calculated by our models. A more detailed chemical network is currently under development. The electron fraction has been computed using a simplified version of the reaction scheme of Umebayashi & Nakano (1990), where a generic molecular ion mH⁺ is formed via proton transfer with H₃⁺, and it is destroyed by dissociative recombination with electrons and recombination on grain surfaces (using rates from Draine & Sutin 1987). Dust grains follow a Mathis et al. (1977; MRN) size distribution, but the minimum size has been increased to  cm to simulate possible dust coagulation, as in the best-fit model of the prestellar core L1544 shown by Vastel et al. 2006 (see also Flower et al. 2005 and Bergin et al. 2006). The initial abundance of metals (assumed to freeze out with a rate similar to that of CO) is 10^-6 (see McKee 1989).

Table 6:   Model results.

The cosmic ray ionization rate is fixed at  s^-1 (van der Tak & van Dishoeck 2000). The adopted CO binding energy is 1100 K, a weighted mean of the CO  binding energy on icy mantles (1180 K; Collings et al. 2003) and CO mantles (885 K; Öberg et al. 2005), assuming that solid water is about four times more abundant than solid CO. The models are computed for envelope lifetimes of 10⁶ yr, although solutions do not appreciably change after 10⁵ yr.

In Figs. 11 to 14, we show the fits for Cep E-mm. We adopt this protostar as a fiducial example because it is the one with the best spatial sampling of the envelope, its geometry looks spherical, and the contribution of the surrounding molecular cloud is negligible. It is impossible to reproduce the observations towards Cep E-mm using the standard chemical model. With the standard model, i.e., assuming a CO binding energy of 1100 K, we are able to reproduce the qualitative behavior of the C¹⁸O and N₂H⁺ emission (see Fig. 11), but the model reproduces the line integrated intensities poorly. The line intensity predictions are a factor of 2-4 higher than the observations (see Fig. 11). Increasing the CO binding energy (e.g. assuming that a large fraction of CO is trapped in water ice) does not improve the fit (see model 2 in Table 6 and Fig. 12).

The abundance profiles given in Table 5 suggest that the fit to the C¹⁸O emission would improve by lowering the CO abundance in the inner part of the core. As pointed out there, this low CO abundance has some physical justification. One possibility is that a significant fraction of CO is converted into CH₃OH on the grain surfaces before evaporation, in agreement with observations of solid methanol along the line of sight toward embedded young stellar objects (e.g. Boogert et al. 2008). We mimic this situation with model 3. In this model, only 10% of the CO is released back to the gas phase at the CO evaporation temperature. Another possibility is that CO is either destroyed by the stellar UV radiation and/or X-rays close to the star or transformed into more complex molecules via hot core chemistry. This case corresponds to our model 4, where only 10% of the CO survives when the temperature is higher than 100 K. The two models, 3 and 4, fit the C¹⁸O and N₂H⁺ emission better than the standard model (see Figs. 13 and 14), with model 3 being a slightly better fit to the Cep E-mm observations. Thus, we applied these most successful models, model 3 and 4, to all the other sources and have obtained reasonable fits (line integrated intensities fitted within a factor of 2) to the C¹⁸O and N₂H⁺ emission in all of them. For both IC 1396 N and Serpens-FIRS 1, model 4 gives a better fit, while model 3 is better for the rest.

While we obtain reasonable fits for C¹⁸O and N₂H⁺, our models do not succeed in reproducing the N₂D⁺ data towards Serpens-FIRS 1 and IC 1396 N. In fact, the chemical models predict intensities higher by a factor of 10 than the observed intensities for these sources. All the considered models predict that the spatial distribution of N₂D⁺ is similar to that of N₂H⁺ and that the deuteration fraction [N₂D⁺]/[N₂H.01 in most of the envelope. This is not true for our low-deuterated sources, Serpens-FIRS 1 and IC 1396 N, in which the [N₂D⁺]/[N₂H⁺] ratio is a few 0.001. This discrepancy could have different origins. First of all, the models assume a spherical geometry. It is clear that bipolar outflows have excavated large cavities in these protostellar envelopes (see e.g. Fuente et al. 2009). The walls of these cavities are warmed by the stellar UV radiation and shocks. Moderate temperatures, UV radiation, and shocks would lower the abundance of N₂H⁺ and change the [N₂D⁺]/[N₂H⁺] ratio. In fact, the large linewidths observed (1.5 km s^-1) are hard to maintain without shocks and dissipation. The role of UV radiation and shocks is expected to be greater in IM Class 0 objects than in low-mass ones.

This simple model ignores many other important parameters that could decrease the deuterium fractionation in warmer sources. These include (i) an increased abundance of atomic O in the gas phase, possibly coming from the release of water from the icy dust mantles. Atomic oxygen is an efficient destruction partner for all the H₃⁺ isotopologues, and lowers the D-fractionation, as discussed by Caselli et al. (2002). (ii) The ortho-to-para H₂ ratio, which in systems out of equilibrium (such as free-falling) can exceed the steady-state value by more than an order of magnitude (Flower et al. 2006; Pagani et al. 2009). A larger fraction of ortho-H₂ leads to a lower D-fractionation, given that the more energetic ortho-H₂ can more easily drive the proton-deuteron exchange reactions (e.g. H₃⁺ + HD H₂D⁺ + H₂) backward and reduce the H₂D⁺/H₃⁺ abundance ratio (Gerlich et al. 2002). (iii) A higher ionization rate, , which may be due to the presence of X-rays. A larger implies a larger electron fraction and thus a higher dissociative recombination rates for molecular ions, including H₃⁺ isotopologues (e.g. Caselli et al. 2008). (iv) The presence of small grains (in particular PAHs; see discussion in Caselli et al. 2008, Sect. 5.1) may arise from the interaction of outflow lobes and UV radiation with the molecular envelope. The associated shocks will partially destroy dust grains along the way (e.g. Jones et al. 1996; Caselli et al. 1997; Guillet et al. 2009) and increase the surface area for recombination of molecular ions onto dust grains.

Because of the low angular resolution of the present observations, it is hard to disentangle the influences of the above parameters. Both more detailed observations and more comprehensive models are needed for a deeper understanding of the chemical and physical evolution of the envelopes surrounding young intermediate-mass stars.

8 Summary and conclusions

We carried out a study of the CO depletion and N₂H⁺ deuteration in a sample of representative IM Class 0 protostars. Our results can be summarized as follows.

• We observed the millimeter lines of C¹⁸O, C¹⁷O, N₂H⁺, and N₂D⁺ using the IRAM 30m telescope in a sample of 7 Class 0 and 2 Class I IM stars. We have found a clear evolutionary trend that differentiates Class 0 from Class I sources. While the emission of the N₂H⁺ 1 0 peaks towards the star position in Class 0 protostars, it surrounds the FIR sources in the case of Class I stars. This occurs because the recently formed star has heated and disrupted the parent core in the case of Class I objects. The deuterium fractionation, R₂, is low, below a few 0.001 in all Class I sources. There is, however, a wide dispersal in the values of R₂ in Class 0 sources ranging from a few 0.001 to a few 0.01. This at least two orders of magnitude greater than the elemental value in the interstellar medium, although a factor of 10-100 lower than in prestellar clumps.

  It is impossible, however, to establish an evolutionary trend among Class 0 sources based on simple parameters such as the average CO depletion and the average N₂H⁺ deuterium fractionation. This stems from the complexity of these regions (multiplicity) and the limited angular resolution of our observations, which prevents us from tracing the inner regions of the envelope. Interferometric observations are required to provide a more precise picture of the evolutionary stage of these objects.

• We used a radiative transfer code to derive the C¹⁸O, N₂H⁺, and N₂D⁺ radial abundance profiles in 5 IM Class 0 stars. In particular, we fit the C¹⁸O 1 0 maps by assuming that the C¹⁸O abundance decreases inwards within the protostellar envelope until the gas and dust reach the CO evaporation temperature, . Our observational data are better fit with values of  K, consistent with the binding energy of 1100 K, which is the measured in the laboratory for a CO-CO matrix.

• We determined the chemistry of the protostellar envelopes using the model by Caselli et al. (2002). A spherical envelope and steady-state chemical model cannot account for the observations. We had to introduce modifications to better fit the C¹⁸O and N₂H⁺ maps. In particular the CO abundance in the inner envelope seems to be lower than the canonical value. This could be due to the conversion of CO into CH₃OH on the grain surfaces, the photodissociation of CO by the stellar UV radiation, or even geometrical effects. Likewise, we have problems fitting the low values of the deuterium fractionation (a few 0.001) measured for some Class 0 IMs. Several explanations have been proposed to account for this discrepancy.

Acknowledgements

This Paper was partially supported by MICINN, within the program CONSOLIDER INGENIO 2010, under grant ``Molecular Astrophysics: The Herschel and ALMA Era - ASTROMOL'' (ref.: CSD2009-00038)

Appendix A:

Below we discuss the details of the modeling for each individual source.

A.1 Serpens-FIRS 1

Serpens-FIRS 1, located near the center of the Serpens main core, is the most luminous object embedded in the cloud. Several continuum studies lead to its classification as a Class 0 source with a bolometric luminosity estimated to range from 46  to 84  (Harvey et al. 1984; Casali et al. 1993; Hurt & Barsony 1996; Larsson et al. 2000). The latest estimate of its luminosity (see Table 4) suggests that Serpens-FIRS 1 is on the low mass/IM borderline. Serpens-FIRS 1 drives a molecular outflow that is oriented at a position angle of 50(Rodríguez et al. 1989). C10 modeled this source as a sphere with an outer radius of 5900 AU, 25.5'' at the distance of Serpens. Convolving this size with the observational beams, one would expect a source diameter of 55'', 57'', and 53'' for the C¹⁸O 1 0, N₂H⁺ 1 0, and N₂D⁺ 2 1 maps, respectively. In our maps, the emission is much more extended (>80''), showing that the dense molecular cloud significantly contributes to the molecular emission (see also Fig. 6). To mimic the molecular cloud component, we increased the outer radius in the Crimier et al. profile. We assumed that the radius of the protostellar envelope is 15 000 AU, and the dust density and temperature smoothly vary between the last point of the Crimier et al. profile (r=5900 AU,  cm^-3,  K) and the values assumed at r=15 000 AU, which are  cm^-3 and  K. This profile was used for our fitting.

The integrated intensity map of C¹⁸O cannot be fit with a constant abundance profile. Thus, we decided to assume a step function for the C¹⁸O abundance: (i) a constant abundance, X₀, that is expected to be close to the canonical value ( ) for radii lower than a given radius, R₀ and (ii) an abundance, X₁, that is expected to be <X₀for larger radii. Thus defined R₀, is the evaporation radius of CO. The values of X₀, R₀, and X₁ were fit by the model. This approximation was still not sufficient to produce a good fit. We had to add another step, and two new variables, R₁ and X₂, in the C¹⁸O abundance profile. The best fit is shown in Table 5. We have (i) a warm region (R<2000 AU) with an C¹⁸O abundance ; (ii) an intermediate layer with a high value for the C¹⁸O depletion, , and (iii) an external layer (R>6000 AU), which is essentially the molecular cloud component, in which the C¹⁸O abundance is close to the canonical value again.

We followed the same procedure for N₂H⁺. In this case the emission extends to the NE and greatly differs from spherical symmetry. For this reason we masked the NE quadrant in our fitting (see Fig. 4). Again, we conclude that the N₂H⁺ abundance has a standard value of a few 10^-10 in the inner region (R< 3000 AU) and in the molecular cloud, but decreases by at least a factor of 10 in the region between them.

Serpens-FIRS 1 is one of the three sources where we were able to fit the N₂D⁺ 2 1 emission. It is clear from the map morphology that the N₂D⁺ emission is dominated by the molecular cloud component. Similar to the case of N₂H⁺, we masked the NE quadrant in our fitting. We found X(N₂D in the molecular cloud, and obtained an upper limit of X(N₂D for the protostellar core. Thus we have a deuterium fractionation of 0.01 in the molecular cloud, and at least an order of magnitude lower in the core (see Fig. 7).

A.2 Cep E

Cep E-mm was cataloged as a Class 0 protostar by Lefloch et al. (1996). Cep E-mm was observed with IRAM 30 m (Lefloch et al. 1996; Chini et al. 2001), SCUBA (Chini et al. 2001), ISO (Froebrich et al. 2003), and Spitzer (Noriega-Crespo et al. 2005). All these studies confirm the Class 0 status of Cep E-mm and constrain the source total mass and bolometric luminosity in the range of 7-25  and 80-120 , respectively. A bipolar molecular outflow, first reported by Fukui et al. (1989), is associated with Cep E-mm. The H₂ and [FeII] study by Eisloffel et al. (1996) shows a quadrupolar outflow morphology suggesting that the driving source is a binary.

C10 modeled this source as a sphere with an outer radius of 35 800 AU, 49'' at the distance of Cepheus (see Table 4). Convolving this size with the observational beams, one would expect a source diameter of 108'', 110'', and 100'' for the C¹⁸O 1 0, N₂H⁺ 1 0, and N₂D⁺ 2 1 maps, respectively. These sizes agree with those found in our maps. Of our sources, Cep E has the largest envelope diameter versus HPBW ratio and thus the best spatial sampling of the varying chemical conditions in its protostellar envelope (see Table 5).

In this source, we fit the emission with a constant and close-to-standard abundance of for radii below 3500 AU, and a power-law variation of the C¹⁸O abundance for larger radii (see Table 5). The high value of , 10, would occur close to R=3500 AU. The same kind of profile was fit for the N₂H⁺ abundance. In this case, the quality of the N₂D⁺ 2 1 line was not good enough to fit the abundance profile. We derived an averaged abundance across the envelope of by fitting the spectrum observed towards the center position.

A.3 IC 1396 N

IC 1396 N is the globule associated with IRAS 21391+5802. It has strong submillimeter and millimeter continuum emission (Wilking et al. 1993; Sugitani et al. 2000; Codella et al. 2001), high density gas (Serabyn et al. 1993; Cesaroni et al. 1999; Codella et al. 2001; Beltrán et al. 2004), and water maser emission (Felli et al. 1992; Tofani et al. 1995; Patel et al. 2000). IC 1396 N is thus an active site of star formation. Using BIMA interferometric millimeter observations, Beltrán et al. (2002) detected three sources (BIMA 1, BIMA 2, and BIMA 3) deeply embedded within the globule. Among the three, BIMA 2 has the strongest millimeter emission and is associated with an energetic E-W bipolar outflow (Codella et al. 2001 Beltrán 2002, 2004). Recent studies by Neri et al. (2007) and Fuente et al. (2007) using the Plateau de Bure interferometer show that BIMA 2 is itself a protocluster composed of 3 cores.

C10 modeled BIMA 2 as a sphere with an outer radius of 29 600 AU, 39'' at the distance of IC 1396 N (see Table 4). Together with Cep E-mm, it is our best-sampled protostellar envelope (large envelope diameter to HPBW ratio). Interferometric observations published by Fuente et al. (2009) reveal that the outflow has already eroded a large biconical cavity in the molecular cloud. This is very likely the reason for the scarce C¹⁸O 1 0 and N₂H⁺ 1 0 emission to the east of the source. We have therefore masked this region to have a more accurate description of the toroidal envelope (see Fig. 4). This masking does not affect our results significantly (more than a factor of 2 in the abundances).

Similarly to Cep-E, the C¹⁸O emission was better fit with a constant and close-to-standard abundance of in the inner region (R<5500 AU) and a power-law variation of the C¹⁸O abundance for larger radii (see Table 5). The highest value of is 5, and it would occur close to R=5500 AU. The N₂H⁺ emission was better fit with a two-step function, for R<10 000 AU, for R> 15 000 AU, and in between. The last step could be due to the vicinity of the sources BIMA 3 and BIMA 2 that heat the outer part of the envelope. These sources are not considered in the (n-T) fit by C10. In IC 1396 N, we have not detected the N₂D⁺ 2 1 line. We derived an upper limit to the N₂D⁺ abundance of < . The lower value of and the non-detection of N₂D⁺ is consistent with this source being a warmer and more evolved object.

A.4 CB3

The CB3 Bok globule is located at 2.5 kpc (Launhardt & Henning 1997; Wang et al. 1995). CB3-mm is the brightest millimeter source of the globule, first detected by Launhardt & Henning (1997) and subsequently observed in the sub-millimeter by Huard et al. (2000). Yun & Clemens (1994) detected a molecular bipolar outflow in CO, elongated in the NE-SW direction, associated with H₂O masers (de Gregorio-Monsalvo et al. 2006). This outflow has been mapped in various molecular lines by Codella & Bachiller (1999), who concluded that it originates from CB3-mm. The same authors concluded that CB3-mm is probably a Class 0 source.

CB3 is different from all the other sources in our sample. In this object, the 24 m sources are not spatially coincident with the 850 m emission peak. CB3-mm hosts two 24 m sources separated by 12'' and neither of them is spatially coincident with the column density peak. The column density peak, better traced by the 850 m continuum emission, is located in between and almost equidistant from the two 24 m sources (see Fig. 4). C10 modeled the SED of this source by adding up the flux of the two 24 m sources. They fit the spatial distribution of the 450 m and 850 m maps as a sphere with an outer radius of 103 000 AU (i.e. 41'') (see Table 4). The radial density power law in this source is steeper, 2, than in the others, with low densities (a few 10³ cm^-3) in the outer envelope, >50 000 AU. This means that, although the protostellar envelope is very large, most of the emission comes from the inner 50 000 AU.

We fit the C¹⁸O emission with a step function. The abundance is constant at for radii less than 25 000 AU, decreases to < ( ), and then increases to for R>60 000 AU (see Table 5). The large C¹⁸O abundance at large radii is very likely caused by the modeled extended low-density envelope providing a poor approximation. More likely, there are dense clumps immersed in a lower density cloud. The N₂H⁺ emission is fit with a constant abundance of . The N₂D⁺ emission is fit with a constant abundance of , implying an average [N₂D⁺]/[N₂H⁺]=0.05.

A.5 NGC 7129-FIRS 2

NGC 7129 is a reflection nebula located in a complex and active star-forming site at a distance of 1250 pc (Hartigan & Lada 1985; Miranda et al. 1993). NGC 7129 FIRS 2 is not detected at optical or near-infrared wavelengths. Its position is spatially coincident with a ¹³CO column density peak (Bechis et al. 1978) and a high-density NH₃ cloudlet (Guesten & Marcaide 1986), and it is close to an H₂O maser (Rodríguez et al. 1980). NGC 7129-FIRS 2 has been classified as a Class 0 IM protostar by Eiroa et al. (1998), who carried out a multi-wavelength study of the continuum emission from 25 to 2000 m. Edwards & Snell (1983) detected a bipolar CO outflow associated with FIRS 2. The interferometric study by Fuente et al. (2001) reveals that this outflow presents a quadrupolar morphology. In fact, the outflow seems to be the superposition of two flows, FIRS 2-out 1 and FIRS 2-out 2, likely associated with FIRS 2 and a more evolved star (FIRS 2-IR), respectively. Fuente et al. (2005a,b) carried out a complete chemical study of FIRS 2 providing the first detection of hot core in an IM Class 0. Based on all these studies, FIRS 2 is considered the youngest IM object known at present.

C10 modeled NGC 7129-FIRS 2 as a sphere with an outer radius of 18 600 AU, 15'' at the distance of NGC 7129 (see Table 4). Since the source is very compact and barely resolved in the maps of the C¹⁸O 1 0 and N₂H⁺ 1 0 lines, we fit only an average abundance in for these species. The average CO depletion factor is , consistent with the depletion factor found by Fuente et al. (2005a) based on H¹³CO⁺ observations. The protostellar envelope is, however, resolved at the frequency of the N₂D⁺ 3 2 line. Since previous interferometric observations published by Fuente et al. (2005a) reveal that the N₂D⁺ emission is absent towards the hot core, we fitted the N₂D⁺ 3 2 emission assuming that it arises in a ring. The best fit is for a ring with inner radius of 9000 AU and an outer radius of 11 000 AU. Within the ring the N₂D⁺ abundance is , and the deuterium fractionation, 0.1, is the typical value found in prestellar cores.

References

Footnotes

... code

The model is called DataCube and a link to install it is available upon request.

... database

LAMDA database is available at http://www.strw.leidenuniv.nl/ moldata/.

All Tables

Table 1:   Selected sample.

Table 2:   Description of the observations.

Table 3:   LTE column densities.

Table 4:   Temperature-density (T-n) profiles from Crimier et al. (2010).

Table 5:   Abundance profiles derived from the observations.

Table 6:   Model results.

All Figures

---

Copyright ESO 2010