# Chemical study of intermediate-mass (IM) Class 0 protostars. CO depletion and N2H+ deuteration

T. Alonso-Albi, A. Fuente, N. Crimier, P. Caselli, C. Ceccarelli, D. Johnstone, P. Planesas, J. R. Rizzo, F. Wyrowski, M. Tafalla, et al.

► To cite this version:

T. Alonso-Albi, A. Fuente, N. Crimier, P. Caselli, C. Ceccarelli, et al., Chemical study of intermediate-mass (IM) Class 0 protostars. CO depletion and N2H+ deuteration. Astronomy & Astrophysics - A&A, 2010, 518, pp.A52. hal-00631908

HAL Id: hal-00631908

https://hal.science/hal-00631908v1

Submitted on 30 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Chemical study of intermediate-mass (IM) Class 0 protostars
## CO depletion and $N_2H^+$ deuteration

T. Alonso-Albi<fnref n="1" />, A. Fuente<fnref n="1" />, N. Crimier<fnref n="2" />, P. Caselli<fnref n="3" />, C. Ceccarelli<fnref n="2" />, D. Johnstone<fnref n="4" />,<fnref n="5" />, P. Planesas<fnref n="1" />,<fnref n="6" />, J. R. Rizzo<fnref n="7" />, F. Wyrowski<fnref n="8" />, M. Tafalla<fnref n="1" />, B. Lefloch<fnref n="2" />, S. Maret<fnref n="2" />, and C. Dominik<fnref n="9" />

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_2H^+$ 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 ($d \le 1$ 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

<footnote n="1">¹ Observatorio Astronómico Nacional (OAN, IGN), Apdo 112, 28803 Alcalá de Henares, Spain
e-mail: t.alonso@oan.es</footnote>

<footnote n="2">² Laboratoire d'Astrophysique, Observatoire de Grenoble, BP 53, 38041 Grenoble Cedex 9, France</footnote>

<footnote n="3">³ School of Physics & Astronomy, E.C. Stoner Building, The University of Leeds, Leeds LS2 9JT, UK</footnote>

<footnote n="4">⁴ Department of Physics & Astronomy, University of Victoria, Victoria, BC, V8P 1A1, Canada</footnote>

<footnote n="5">⁵ National Research Council of Canada, Herzberg Institute of Astrophysics, 5071 West Saanich Road, Victoria, BC, V9E 2E7, Canada</footnote>

<footnote n="6">⁶ Joint ALMA Observatory, El Golf 40, Las Condes, Santiago, Chile</footnote>

<footnote n="7">⁷ Centro de Astrobiología (CSIC/INTA), Laboratory of Molecular Astrophysics, Ctra. Ajalvir km. 4, 28850, Torrejón de Ardoz, Spain</footnote>

<footnote n="8">⁸ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>

<footnote n="9">⁹ Anton Pannekoek Astronomical Institute, University of Amsterdam, PO Box 94249, 1090 GE Amsterdam, The Netherlands</footnote>
<table><caption>**Table 1.** Selected sample.</caption><thead><tr><th>Object</th><th>RA(2000)</th><th>Dec(2000)</th><th>Lum. (L<sub>☉</sub>)</th><th>d</th><th></th></tr></thead><tbody><tr><td>Serpens-FIRS 1</td><td>18:29:49.6</td><td>+01:15:20.6</td><td>33</td><td>230</td><td>Maps in C<sup>17</sup>O, C<sup>18</sup>O, N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup></td></tr><tr><td>Cep E-mm</td><td>23:03:13.1</td><td>+61:42:26.0</td><td>100</td><td>730</td><td>Maps in C<sup>17</sup>O, C<sup>18</sup>O, N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup></td></tr><tr><td>L1641 S3 MMS 1</td><td>05:39:55.9</td><td>-07:30:28.0</td><td>67</td><td>500</td><td>Strip in C<sup>17</sup>O and C<sup>18</sup>O, maps in N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup></td></tr><tr><td>IC1396N</td><td>21:40:41.7</td><td>+58:16:12.8</td><td>150</td><td>750</td><td>Maps in C<sup>17</sup>O, C<sup>18</sup>O, N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup></td></tr><tr><td>CB3</td><td>00:28:42.7</td><td>+56:42:07.1</td><td>1000</td><td>2500</td><td>Maps in C<sup>17</sup>O, C<sup>18</sup>O, N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup></td></tr><tr><td>OMC2-FIR4</td><td>05:35:26.7</td><td>-05:10:00.5</td><td>1000</td><td>450</td><td>Maps in C<sup>17</sup>O, C<sup>18</sup>O, N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup></td></tr><tr><td>NGC 7129-FIRS 2</td><td>21:43:01.7</td><td>+66:03:23.6</td><td>500</td><td>1250</td><td>Strip in C<sup>17</sup>O, maps in C<sup>18</sup>O, N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup></td></tr><tr><td>S140</td><td>22:19:18.1</td><td>+63:18:54.6</td><td>10<sup>4</sup></td><td>910</td><td>N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> maps</td></tr><tr><td>LkHα 234</td><td>21:43:06.8</td><td>+66:06:54.4</td><td>500</td><td>1250</td><td>Maps in C<sup>18</sup>O, N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup></td></tr></tbody></table>

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<sup>18</sup>O, C<sup>17</sup>O, N<sub>2</sub>H<sup>+</sup>, and N<sub>2</sub>D<sup>+</sup> in the same sample. Our goal is to investigate the CO depletion and N<sub>2</sub>H<sup>+</sup> deuteration in these Class 0 YSOs. For comparison, we also include 2 Class I objects, LkHα 234, and S140.

## 2. Observational strategy

Our selection was made to have a representative sample of Class 0 IM YSOs, including targets with different luminosities (40–10<sup>3</sup> L<sub>☉</sub>) 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" × 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<sup>17</sup>O and C<sup>18</sup>O. We did not observe C<sup>17</sup>O and C<sup>18</sup>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<sub>2</sub>H<sup>+</sup> *J* = 4 → 3 line (freq = 372.6725 GHz) were carried out towards Cep E-mm, IC 1396 N,

<table><caption>**Table 2.** Description of the observations.</caption><thead><tr><th>Line</th><th>Freq. (MHz)</th><th>HPBW</th><th>η<sub>MB</sub></th><th>Tel</th></tr></thead><tbody><tr><td>C<sup>18</sup>O 1→0</td><td>109 782.2</td><td>22"</td><td>0.75</td><td>IRAM</td></tr><tr><td>C<sup>18</sup>O 2→1</td><td>219 560.3</td><td>11"</td><td>0.55</td><td>IRAM</td></tr><tr><td>C<sup>17</sup>O 1→0</td><td>112 359.3</td><td>22"</td><td>0.74</td><td>IRAM</td></tr><tr><td>C<sup>17</sup>O 2→1</td><td>224 714.4</td><td>11"</td><td>0.54</td><td>IRAM</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup> 1→0</td><td>93 173.2</td><td>26"</td><td>0.77</td><td>IRAM</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup> 4→3</td><td>372 672.5</td><td>13"</td><td>0.63</td><td>JCMT</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup> 2→1</td><td>154 217.1</td><td>16"</td><td>0.67</td><td>IRAM</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup> 3→2</td><td>231 321.7</td><td>11"</td><td>0.52</td><td>IRAM</td></tr></tbody></table>

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" × 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 *J* = 4 → 3 line are shown in Fig. 3.

In this paper, we also use the C<sup>18</sup>O, N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> data towards NGC 7129-FIRS 2 and LkHα 234, which has already been published by Fuente et al. (2005b).

## 3. Results

The spectra of the C<sup>18</sup>O 1→0, C<sup>17</sup>O 1→0, N<sub>2</sub>H<sup>+</sup> 1→0 and N<sub>2</sub>D<sup>+</sup> 2→1 lines towards the star position are shown in Figs. 1 and 2. The integrated intensity maps of the C<sup>18</sup>O 1→0, N<sub>2</sub>H<sup>+</sup> 1→0, and N<sub>2</sub>D<sup>+</sup> 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<sub>2</sub>H<sup>+</sup> 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<sub>2</sub>H<sup>+</sup> emission are offset from the far-IR source. For the Herbig Be star LkHα 234, the emission peak of the N<sub>2</sub>H<sup>+</sup> 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<sub>2</sub>H<sup>+</sup> emission peaks in an arc-shaped feature surrounding the sources IRS 1, 2 and 3. In Class I sources, either the N<sub>2</sub>H<sup>+</sup> 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<sup>18</sup>O 1→0 line usually presents a different morphology from that of N<sub>2</sub>H<sup>+</sup> 1→0 line. In CB3 and IC 1396N, the emission in the C<sup>18</sup>O line has an elongated shape, much more extended than that of N<sub>2</sub>H<sup>+</sup> (see Fig. 4). In OMC2 FIR 4 and NGC 7129-FIRS 2, the emission of the C<sup>18</sup>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<sup>17</sup>O and C<sup>18</sup>O lines peak at the star position and present a morphology similar to that of the N<sub>2</sub>H<sup>+</sup> emission.
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Examples of the single-pointing observations towards the Serpens-FIRS 1, Cep E-mm, L1641 S3 MMS1, and IC 1396 N. All the spectra were observed with the 30 m telescope. The intensity scale is the main brightness temperature.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** The same as Fig. 1 for CB3, OMC2 FIR 4, NGC 7129-FIRS 2, and S140.</figcaption></figure>

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
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Spectra of the N<sub>2</sub>H<sup>+</sup> 4→3 line towards NGC 7129–FIRS 2, Cep E-mm, IC 1396 N, and L1641 S3 MMS1.</figcaption></figure>

the case of a very young object or to photodissociation for a borderline Class 0/I. The radial profiles of the N<sub>2</sub>D<sup>+</sup> 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<sup>18</sup>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<sup>18</sup>O and C<sup>17</sup>O are very similar, and we have better signal-to-noise spectra for the more intense C<sup>18</sup>O lines. For this reason, we derived the C<sup>17</sup>O column density by assuming optically thin emission and the rotation temperature derived from the C<sup>18</sup>O data. The molecule N<sub>2</sub>H<sup>+</sup> presents hyperfine splitting. This allows us to estimate the line opacity directly from the hyperfine line ratios. We derived the total N<sub>2</sub>H<sup>+</sup> column density from the opacity of the N<sub>2</sub>H<sup>+</sup> 1→0 line, and assumed the rotation temperatures derived from N<sub>2</sub>D<sup>+</sup> data when there was no N<sub>2</sub>H<sup>+</sup> 4→3 observations. The same rotation temperature was always used for both N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup>. 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<sup>18</sup>O, C<sup>17</sup>O, N<sub>2</sub>H<sup>+</sup>, and N<sub>2</sub>D<sup>+</sup> column densities. We also show the N(C<sup>18</sup>O)/N(C<sup>17</sup>O) (hereafter R<sub>1</sub>) and the N(N<sub>2</sub>D<sup>+</sup>)/N(N<sub>2</sub>H<sup>+</sup>) (hereafter R<sub>2</sub>) ratios. For the whole sample, we obtain a value for R<sub>1</sub> around 3.3+/-0.3, the expected value for optically thin emission. In the worst case, R<sub>1</sub> = 2.5, the correction to the C<sup>18</sup>O column density because of the opacity is only a factor ~1.3. Thus, opacity effects are not important in our C<sup>18</sup>O column density estimates.

**4.1. N<sub>2</sub>D<sup>+</sup> and N<sub>2</sub>H<sup>+</sup>**

The deuterium fractionation of N<sub>2</sub>H<sup>+</sup> (R<sub>2</sub>) 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<sub>2</sub> ranging from 0.005 to 0.014. These values are 3 orders of magnitude higher than 10<sup>-5</sup>, the elemental value in the interstellar medium

(Oliveira et al. 2003). Nevertheless, the R<sub>2</sub> values are a factor of 10 lower than those found in prestellar clumps by Crapsi et al. (2005). According to the values of R<sub>2</sub>, we can classify our sources into two groups: (i) highly deuterated sources that have values of R<sub>2</sub> > 0.01 and (ii) moderately deuterated sources with R<sub>2</sub> < 0.01. Cep E-mm, CB3, and NGC 7129–FIRS 2 belong to the first group. Assuming that the N(N<sub>2</sub>D<sup>+</sup>)/N(N<sub>2</sub>H<sup>+</sup>) 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<sub>2</sub>D<sup>+</sup> 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<sup>18</sup>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<sub>2</sub>H<sup>+</sup> deuterium fractionation as a function of the N(C<sup>18</sup>O)/N(N<sub>2</sub>H<sup>+</sup>) 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<sup>18</sup>O)/N(N<sub>2</sub>H<sup>+</sup>) 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.

**5. Radiative transfer model**

We utilized a general ray-tracing radiative transfer code<fnref n="1" /> to derive the fractional abundance profiles of C<sup>18</sup>O, N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> 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 32 × 32 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<sup>-1</sup>, 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

<footnote n="1">¹ The model is called DataCube and a link to install it is available upon request.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Continuum maps at 850 µm from C10 (left column), and integrated intensity maps of the C¹⁸O 1→0, N₂H⁺ 1→0 and N₂D⁺ 2→1 (3→2 in the case of NGC 7129-FIRS 2) lines towards the observed Class 0 IMs. In each panel, we have drawn the area and the number of rings considered in our model fitting. The emission peak at several wavelengths is shown with different marks. Contour levels start at 3σ level and are **a)** 0.5, 1.6, 2.7, 3.8, 4.9, and 5.4 Jy/beam; **b)** 1 to 9 K km s⁻¹ by steps of 1 K km s⁻¹; **c)** 3 to 9 K km s⁻¹ by steps of 3 K km s⁻¹; **d)** 0.5 to 1.5 K km s⁻¹ by steps of 0.5 K km s⁻¹; **e)** 0.16, 0.49, 0.81, 1.14, and 1.46 Jy/beam; **f)** 0.5 to 3 K km s⁻¹ by steps of 0.5 K km s⁻¹; **g)** 0.5 to 2.5 K km s⁻¹ by steps of 0.5 K km s⁻¹; **h)** 0.3, 0.9, 1.5, 2.1, and 2.8 Jy/beam; **i)** 2.0 to 8.0 K km s⁻¹ by steps of 2.0 K km s⁻¹; **j)** 3.0 to 9.0 K km s⁻¹ by steps of 3.0 K km s⁻¹; **k)** 0.07, 0.20, 0.33, 0.46, and 0.59 Jy/beam; **l)** 2.0, 4.0 K km s⁻¹; **m)** 1.0 to 3.0 K km s⁻¹ by steps of 1.0 K km s⁻¹; **n)** 1.0, 2.5, 4.0, 5.5, and 7.0 Jy/beam; **o)** 0.5 to 3.5 K km s⁻¹ by steps of 0.5 K km s⁻¹; **p)** 4.0 to 24.0 K km s⁻¹ by steps of 4.0 K km s⁻¹; **q)** 2.1, 6.3, 10.6, 14.8, and 19.0 Jy/beam; **r)** 2.0 K km s⁻¹; **s)** 3.0, 6.0 K km s⁻¹; **t)** 0.5 K km s⁻¹.</figcaption></figure>
A&A 518, A52 (2010)

<figure><img src="image_6.png" /></figure>

**Fig. 5.** Integrated intensity maps of the C<sup>18</sup>O 1→0 and N<sub>2</sub>H<sup>+</sup> 1→0 lines towards LkHα 234 and S140. Contour levels are: **a)** 0.07, 0.15, 0.23, 0.31, 0.39, 0.46, 0.54, 0.62, and 0.69 Jy/beam; **b)** 2 to 5 K km s<sup>-1</sup> by steps of 1 K km s<sup>-1</sup>; **c)** 1 to 5 K km s<sup>-1</sup> in steps of 1 K km s<sup>-1</sup>; **d)** 1 to 8.0 Jy/beam in steps of 1 Jy/beam; **e)** 3 to 21 K km s<sup>-1</sup> in steps of 3 K km s<sup>-1</sup>.

<figure><img src="image_7.png" /></figure>

**Fig. 6.** Radially integrated intensity profiles of the C<sup>18</sup>O 1→0, N<sub>2</sub>H<sup>+</sup> 1→0, and N<sub>2</sub>D<sup>+</sup> 2→1 (3→2 for NGC 7129-FIRS 2) lines in our sample of Class 0 IMs. Each point represents the mean value of the line integrated emission in a ring of a given radius. The errors are the rms of the line integrated emission in the ring. The step between two consecutive points, which is equal to the size of each ring, is equivalent to 1/4 of the beam at each frequency.
<table><caption>**Table 3.** LTE column densities.</caption><thead><tr><th>Source</th><th>T<sub>r</sub><br>(K)</th><th>N(C<sup>18</sup>O)<br>(cm<sup>-2</sup>)</th><th>N(C<sup>17</sup>O)<br>(cm<sup>-2</sup>)</th><th>R<sub>1</sub><sup>*</sup></th><th>T<sub>r</sub> (K)<br>(K)</th><th>N(N<sub>2</sub>H<sup>+</sup>)<br>(cm<sup>-2</sup>)</th><th>N(N<sub>2</sub>D<sup>+</sup>)<br>(cm<sup>-2</sup>)</th><th>R<sub>2</sub><sup>*</sup></th></tr></thead><tbody><tr><td>Serpens-FIRS 1</td><td>11</td><td>8.2 &times; 10<sup>15</sup></td><td>3.2 &times; 10<sup>15</sup></td><td>2.6</td><td>11<sup>a</sup></td><td>7.7 &times; 10<sup>13</sup></td><td>3.4 &times; 10<sup>11</sup></td><td>0.004</td></tr><tr><td>Cep E-mm</td><td>15</td><td>2.5 &times; 10<sup>15</sup></td><td>7.6 &times; 10<sup>14</sup></td><td>3.3</td><td>9<sup>b</sup></td><td>1.4 &times; 10<sup>13</sup></td><td>5.8 &times; 10<sup>11</sup></td><td>0.04</td></tr><tr><td>L1641 S3 MMS 1</td><td>10</td><td>3.2 &times; 10<sup>15</sup></td><td>1.3 &times; 10<sup>15</sup></td><td>2.5</td><td>8<sup>b</sup></td><td>3.6 &times; 10<sup>13</sup></td><td>&lt;2.1 &times; 10<sup>11</sup></td><td>&lt;0.006</td></tr><tr><td>IC1396N</td><td>25</td><td>1.1 &times; 10<sup>16</sup></td><td>3.7 &times; 10<sup>15</sup></td><td>3.0</td><td>10<sup>b</sup></td><td>4.1 &times; 10<sup>13</sup></td><td>&lt;4.1 &times; 10<sup>11</sup></td><td>&lt;0.01</td></tr><tr><td>CB3</td><td>14</td><td>5.2 &times; 10<sup>15</sup></td><td>1.6 &times; 10<sup>15</sup></td><td>3.3</td><td>6<sup>a</sup></td><td>1.2 &times; 10<sup>13</sup></td><td>1.0 &times; 10<sup>12</sup></td><td>0.08</td></tr><tr><td>OMC2-FIR4</td><td>22</td><td>4.8 &times; 10<sup>15</sup></td><td>1.2 &times; 10<sup>15</sup></td><td>3.4</td><td>22<sup>a</sup></td><td>2.0 &times; 10<sup>14</sup></td><td>8.7 &times; 10<sup>11</sup></td><td>0.004</td></tr><tr><td>NGC 7129-FIRS 2</td><td>23</td><td>1.9 &times; 10<sup>15</sup></td><td>6.9 &times; 10<sup>14</sup></td><td>2.8</td><td>13<sup>c</sup></td><td>3.8 &times; 10<sup>13</sup></td><td>5.4 &times; 10<sup>11</sup></td><td>0.014</td></tr><tr><td>S140</td><td>50<sup>c</sup></td><td>1.0 &times; 10<sup>16</sup></td><td>2.7 &times; 10<sup>15</sup></td><td>3.7</td><td>50<sup>d</sup></td><td>1.2 &times; 10<sup>14</sup></td><td>&lt;4.5 &times; 10<sup>11</sup></td><td>&lt;0.004</td></tr><tr><td>LkHα 234<sup>b</sup></td><td></td><td></td><td></td><td></td><td>22</td><td>1.0 &times; 10<sup>13</sup></td><td>&lt;2.3 &times; 10<sup>11</sup></td><td>&lt;0.02</td></tr></tbody></table>

**Notes.** (*) $R_1 = N(C^{18}O)/N(C^{17}O)$, $R_2 = N(N_2D^+)/N(N_2H^+)$.

(a) Rotation temperature derived from our N₂D⁺ observations.

(b) Rotation temperature derived from our N₂H⁺ observations.

(c) From Fuente et al. (2005a).

(d) From Minchin et al. (1995).

<table><caption>**Table 4.** Temperature-density ($T-n$) profiles from Crimier et al. (2010).</caption><thead><tr><th>Source</th><th>CB3<br>mm</th><th>Cep- E<br>mm</th><th>IC 1396 N<br>BIMA2</th><th>NGC 7129<br>FIRS 2</th><th>Serpens<br>FIRS 1</th></tr></thead><tbody><tr><td>RA(J2000)</td><td>00:28:42.1</td><td>23:03:12.7</td><td>21:40:41.8</td><td>21:43:01.5</td><td>18:29:49.8</td></tr><tr><td>Dec(J2000)</td><td>56:41:59.4</td><td>61:42:27.4</td><td>58:16:13.5</td><td>66:03:25.0</td><td>01:15:18.4</td></tr><tr><td>Dust optical depth at 100 µm, τ<sub>100</sub></td><td>5.8</td><td>5.0</td><td>1.4</td><td>2.3</td><td>3.0</td></tr><tr><td>Density power law index, α</td><td>2.2</td><td>1.9</td><td>1.2</td><td>1.4</td><td>1.5</td></tr><tr><td>Envelope thickness, r<sub>out</sub>/r<sub>i</sub></td><td>400</td><td>500</td><td>630</td><td>180</td><td>200</td></tr><tr><td>Inner envelope radius, r<sub>in</sub>, (AU)</td><td>260</td><td>70</td><td>50</td><td>100</td><td>30</td></tr><tr><td>Outer envelope radius, r<sub>out</sub>, (AU)</td><td>103 000</td><td>35 800</td><td>29 600</td><td>18 600</td><td>5900</td></tr><tr><td>Radius at T<sub>dust</sub> = 100 K, r<sub>100 K</sub>, (AU)</td><td>700</td><td>223</td><td>180</td><td>373</td><td>102</td></tr><tr><td>H<sub>2</sub> density at r<sub>100 K</sub>, n<sub>0</sub>, (cm<sup>-3</sup>)</td><td>7.5 &times; 10<sup>7</sup></td><td>2.0 &times; 10<sup>8</sup></td><td>4.3 &times; 10<sup>7</sup></td><td>4.4 &times; 10<sup>7</sup></td><td>2.2 &times; 10<sup>8</sup></td></tr><tr><td>Envelope mass, M<sub>env</sub>, (M<sub>☉</sub>)</td><td>120</td><td>35</td><td>90</td><td>50</td><td>5.0</td></tr></tbody></table>

the collisional rates provided by the LAMDA database<fnref n="2" /> (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 ($R < 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 each 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.

<footnote n="2">² LAMDA database is available at http://www.strw.leidenuniv.nl/~moldata/.</footnote>
<table><caption>**Table 5.** Abundance profiles derived from the observations.</caption><thead><tr><th></th><th>X<sub>0</sub></th><th>R<sub>0</sub> (AU)</th><th>X<sub>1</sub></th><th>R<sub>1</sub> (AU)</th><th>X<sub>2</sub></th><th>rms (K km s<sup>-1</sup>)</th><th>Diameter("")/HPBW("")</th></tr></thead><tbody><tr><td><strong>Serpens-FIRS 1</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>C<sup>18</sup>O</td><td>1.4 &times; 10<sup>-7</sup></td><td>2000</td><td>&lt;7 &times; 10<sup>-9</sup></td><td>6000</td><td>1.1 &times; 10<sup>-7</sup></td><td>0.27</td><td>2.3</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup></td><td>4.0 &times; 10<sup>-10</sup></td><td>3000</td><td>&lt;3 &times; 10<sup>-11</sup></td><td>6000</td><td>2.5 &times; 10<sup>-10</sup></td><td>0.29</td><td>1.9</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup></td><td></td><td>&lt;2 &times; 10<sup>-12</sup></td><td></td><td>6000</td><td>1.9 &times; 10<sup>-11</sup></td><td>0.6</td><td>3.1</td></tr><tr><td>T<sub>ev</sub> (K)</td><td>20</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>n<sub>ev</sub><sup>*</sup> (cm<sup>-3</sup>)</td><td>2.4 &times; 10<sup>6</sup></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>Cep E-mm</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>C<sup>18</sup>O</td><td>6.0 &times; 10<sup>-8</sup></td><td>3500</td><td></td><td>1.1 &times; 10<sup>-7</sup> &times; (r/35 800)</td><td></td><td>0.24</td><td>4.4</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup></td><td>2.4 &times; 10<sup>-10</sup></td><td>6000</td><td></td><td>2.5 &times; 10<sup>-10</sup> &times; (r/35 800)</td><td></td><td>0.28</td><td>3.7</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup></td><td></td><td></td><td>7.0 &times; 10<sup>-12</sup></td><td></td><td></td><td></td><td>6.1</td></tr><tr><td>T<sub>ev</sub> (K)</td><td>20</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>n<sub>ev</sub><sup>*</sup> (cm<sup>-3</sup>)</td><td>1.0 &times; 10<sup>6</sup></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>IC 1396 N</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>C<sup>18</sup>O</td><td>6.0 &times; 10<sup>-8</sup></td><td>5500</td><td></td><td>1.25 &times; 10<sup>-7</sup> &times; (r/29 600)</td><td></td><td>0.24</td><td>3.5</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup></td><td>4.2 &times; 10<sup>-10</sup></td><td>10 000</td><td></td><td>1.0 &times; 10<sup>-11</sup></td><td>15 000</td><td>1.8 &times; 10<sup>-10</sup></td><td>3.0</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup></td><td></td><td></td><td>&lt;7 &times; 10<sup>-12</sup></td><td></td><td></td><td></td><td>4.9</td></tr><tr><td>T<sub>ev</sub> (K)</td><td>19</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>n<sub>ev</sub><sup>*</sup> (cm<sup>-3</sup>)</td><td>7.0 &times; 10<sup>5</sup></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>CB3</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>C<sup>18</sup>O</td><td>1.3 &times; 10<sup>-7</sup></td><td>25 000</td><td>&lt;3.0 &times; 10<sup>-8</sup></td><td>60 000</td><td>9.0 &times; 10<sup>-7</sup></td><td>0.24</td><td>3.7</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup></td><td></td><td></td><td>6.5 &times; 10<sup>-10</sup></td><td></td><td></td><td>0.45</td><td>3.1</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup></td><td></td><td></td><td>3.0 &times; 10<sup>-11</sup></td><td></td><td></td><td></td><td>5.1</td></tr><tr><td>T<sub>ev</sub> (K)</td><td>15</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>n<sub>ev</sub><sup>*</sup> (cm<sup>-3</sup>)</td><td>3.0 &times; 10<sup>4</sup></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>NGC 7129-FIRS 2</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>C<sup>18</sup>O</td><td></td><td></td><td>4.0 &times; 10<sup>-8</sup></td><td></td><td></td><td>0.4</td><td>1.4</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup></td><td></td><td></td><td>4.7 &times; 10<sup>-10</sup></td><td></td><td></td><td>0.4</td><td>1.1</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup></td><td>&lt;3.0 &times; 10<sup>-12</sup></td><td>9000</td><td>4.5 &times; 10<sup>-11</sup></td><td>11 000</td><td>&lt;3.0 &times; 10<sup>-12</sup></td><td>0.13</td><td>2.7</td></tr></tbody></table>

**Notes.** (*) Temperature and density at the CO evaporation radius.

Another important source of uncertainty is the degeneracy of the solutions. For instance, in the case of C<sup>18</sup>O we can obtain similar results by varying the molecular abundance in the inner region (X<sub>0</sub>) or the CO evaporation radius (R<sub>0</sub>), 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<sup>18</sup>O abundance profile is X = X<sub>0</sub> for radii less than the evaporation radius, R<sub>0</sub>, and X = X<sub>0</sub> &times; (R/R<sub>out</sub>)<sup>&alpha;</sup> for radii larger than R<sub>0</sub>. For each value of X<sub>0</sub>, we fit R<sub>0</sub> and &alpha;. X<sub>0</sub> is the un-depleted C<sup>18</sup>O abundance, R<sub>0</sub> the evaporation radius, and &alpha; measures the gradient in the C<sup>18</sup>O abundance due to depletion. The canonical abundance of C<sup>18</sup>O in principle depends on the Galactocentric distance (DGC). Following Wilson & Matteucci (1992), the CO abundance is given by

$$X_{(CO)} = 9.5 \times 10^5 \exp(1.105 - (0.13 \times \text{DGC(kpc)})). \quad (1)$$

Following Wilson & Rood (1994), the oxygen isotope ratio <sup>16</sup>O/<sup>18</sup>O depends on DGC according to the relationship <sup>16</sup>O/<sup>18</sup>O = 58.8 &times; DGC(kpc) + 37.1. At the distance of the Sun, 8.5 kpc, X<sub>C<sup>18</sup>O</sub> = 1.7 &times; 10<sup>-7</sup>. This is an average value,

so there may be local effects. In our models we varied X<sub>0</sub> from 1 &times; 10<sup>-7</sup> to 4 &times; 10<sup>-7</sup>, the range of X<sub>0</sub> values that we consider reasonable. For each value of X<sub>0</sub>, we varied R<sub>0</sub> and &alpha; to find the best fit. In Fig. 8 we show the results of our models for Cep E. The best fit is X<sub>0</sub> = 1 &times; 10<sup>-7</sup>, R<sub>0</sub> = 2510 AU, &alpha; = 1 with rms = 0.37 K km s<sup>-1</sup>. 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 X<sub>0</sub> = 1 &times; 10<sup>-7</sup>, the best fit is obtained with &alpha; ~ 1 and R<sub>0</sub> ~ 2100–3100 AU. Assuming X<sub>0</sub> = 2 &times; 10<sup>-7</sup>, we have reasonably good solutions with &alpha; ~ 2 and R<sub>0</sub> ~ 2200–3200 AU. We have not found any acceptable solution for higher values of X<sub>0</sub>. The evaporation radius of C<sup>18</sup>O is thus 2600 &plusmn; 600 AU. These radii correspond to dust temperatures of 20–25 K. Interestingly, we would have a better fit to the observations (rms = 0.24 K km s<sup>-1</sup>) if we allowed the C<sup>18</sup>O abundance in the inner region to fall below 1 &times; 10<sup>-7</sup>. 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<sup>18</sup>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<sub>3</sub>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.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Deuterium fractionation of $N_2H^+$($R_2$) vs. the $N(^{18}O)/N(^{14}N)$ ($N_2H^+$) ratio in the studied YSOs. Inverted empty triangles are upper limits for L1641 S3 MMS1, IC 1396 N, and S140, where we did not detect the $N_2D^+$ 2→1 line.</figcaption></figure>

In Fig. 9 we show the results for IC 1396 N. The best fit is $X_0 = 1 \times 10^{-7}$, $R_0 = 2000$ AU, and $\alpha = 0.5$ with rms = 0.36 K km s⁻¹. We consider that the models with values of rms < 0.5 are acceptable. Following this criterion, we only find solutions for $X_0 = 1 \times 10^{-7}$ and $R_0 < 3000$ AU; i.e., the CO evaporation temperature should be >23 K. Again, the fit can be significantly improved if values of $X_0$ lower than $1 \times 10^{-7}$ 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 $X_0 = 4 \times 10^{-7}$, $R_0 = 4000$ AU, and $\alpha = 0.76$ with rms = 2.0 K km s⁻¹, 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 $T_{ev} \sim 15$ 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 $^{18}$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_2H^+$ and $N_2D^+$ emission profiles have been done by hand. For $N_2H^+$ we tried two options: (i) constant $N_2H^+$ abundance and (ii) an abundance profile with a radial variation similar to that of $^{18}$O. Option (ii) follows the expectation that the $N_2H^+$ 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_2H^+$ observations by assuming a constant abundance. For the rest of sources, the $N_2H^+$ observations are better fit by assuming that the $N_2H^+$ is depleted in the cold regions similarly to CO.

One expects an annular distribution for the $N_2D^+$ abundance, with the maximum $N_2D^+$ abundance in the region with the greatest CO depletion. The morphology observed in the integrated intensity map of the $N_2D^+$ 2→1 line towards Serpens-FIRS 1 suggests that the $N_2D^+$ 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

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Plots of the rms, defined as $\Sigma(I_{model} - I_{obs})^2$ for the grid of models run to reproduce the $^{18}$O emission in Cep E-mm. The $^{18}$O abundance is assumed to be $X = X_0$ for radii less than $R_0$, and $X = X_0 \times (R/R_{out})^\alpha$ for radii larger than $R_0$. The values of $X_0$ are $1.0 \times 10^{-7}$ (panel a)), $2.0 \times 10^{-7}$ (panel b)), $3.0 \times 10^{-7}$ (panel c)), and $4.0 \times 10^{-7}$ (panel d)). The solution with the lowest rms, 0.37, is indicated by a cross. We consider acceptable solutions those with rms less than 0.5. Contour levels are 0.5, 1, 2, and 5 K km s⁻¹.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** The same as Fig. 8 for IC 1396 N. The solution with the lowest rms, 0.36, is indicated by a cross. We consider acceptable solutions those with rms less than 0.5. Contour levels are 0.5, 1, 2, and 5 K km s⁻¹.</figcaption></figure>

clump of $N_2D^+$ 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_2D^+$ 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_2$ 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
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** The same as Fig. 8 for CB 3. The solution with the lowest rms, 2.5, is indicated by a cross. Contour levels are 2.5, 4, 6, and 10 K km s⁻¹.</figcaption></figure>

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 9.5 × 10⁻⁵ (Frerking et al. 1982) and 2 × 10⁻⁵, 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₂) ≤ 2 × 10⁻⁶ (Hily-Blant et al. 2010), and X(N₂) ~ 10⁻⁶ (Maret et al. 2006). Our adopted value is consistent with 25% of the total abundance of N (n(N)/n(H) = 7.9 × 10⁻⁵, 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 5 × 10⁻⁶ 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⁻⁶ (see McKee 1989).

The cosmic ray ionization rate is fixed at $\zeta = 3 \times 10^{-17} \text{ 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 $\approx 10^5$ 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 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
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Left: results of our chemical model for Cep E-mm assuming a binding energy of 1100 K (standard value, model 1). The C<sup>18</sup>O abundance is shown in black, N<sub>2</sub>H<sup>+</sup> in red, and N<sub>2</sub>D<sup>+</sup> in blue. Center: comparison between the predicted C<sup>18</sup>O 1→0 line intensities (continuous line) and the observed values (filled squares). Right: same for N<sub>2</sub>H<sup>+</sup> 1→0.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** The same as Fig. 11 for model 2.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** The same as Fig. 11 for model 3.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** The same as Fig. 11 for model 4.</figcaption></figure>
<table><caption>**Table 6.** Model results.</caption><thead><tr><th>Source</th><th>Model</th><th>CO Binding energy</th><th>χ<sub>C<sup>18</sup>O</sub><sup>a</sup></th><th>χ<sub>N<sub>2</sub>H<sup>+</sup></sub></th><th>χ<sub>N<sub>2</sub>D<sup>+</sup></sub></th><th></th></tr></thead><tbody><tr><td rowspan="2">Serpens-FIRS 1</td><td>Model 3</td><td>1100</td><td>4.37</td><td>3.89</td><td>4.52</td><td>only 10% of CO evaporated</td></tr><tr><td><strong>Model 4</strong></td><td><strong>1100</strong></td><td><strong>2.55</strong></td><td><strong>3.77</strong></td><td><strong>4.46</strong></td><td><strong>only 10% survives for T > 100 K</strong></td></tr><tr><td rowspan="4">Cep E-mm</td><td>Model 1</td><td>1100</td><td>1.48</td><td>4.12</td><td>0.71</td><td></td></tr><tr><td>Model 2</td><td>5000</td><td>1.51</td><td>0.86</td><td>0.45</td><td></td></tr><tr><td><strong>Model 3</strong></td><td><strong>1100</strong></td><td><strong>0.91</strong></td><td><strong>2.26</strong></td><td><strong>1.19</strong></td><td><strong>only 10% of CO evaporated</strong></td></tr><tr><td>Model 4</td><td>1100</td><td>1.42</td><td>2.19</td><td>1.12</td><td>only 10% survives for T > 100 K</td></tr><tr><td rowspan="2">IC 1396 N</td><td>Model 3</td><td>1100</td><td>2.95</td><td>4.45</td><td>>4.3</td><td>only 10% of CO evaporated</td></tr><tr><td><strong>Model 4</strong></td><td><strong>1100</strong></td><td><strong>1.79</strong></td><td><strong>4.35</strong></td><td><strong>>4.0</strong></td><td><strong>only 10% survives for T > 100 K</strong></td></tr><tr><td rowspan="2">CB3</td><td>Model 3</td><td>1100</td><td>1.63</td><td>0.76</td><td>0.31</td><td><strong>only 10% of CO evaporated</strong></td></tr><tr><td>Model 4</td><td>1100</td><td>1.24</td><td>0.86</td><td>0.51</td><td>only 10% survives for T > 100 K</td></tr><tr><td rowspan="2">NGC 7129-FIRS 2</td><td>Model 3</td><td>1100</td><td>0.97</td><td>1.29</td><td>2.05</td><td><strong>only 10% of CO evaporated</strong></td></tr><tr><td>Model 4</td><td>1100</td><td>1.39</td><td>1.48</td><td>1.73</td><td>only 10% survives for T > 100 K</td></tr></tbody></table>

Notes. (a) $\chi$ is defined as $\chi = \sqrt{[\Sigma(I_{model} - I_{obs})^2]/N}$. Note that this is an absolute error and the comparison among values in different sources is not straightforward. See Figs. 11 to 18.

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<sup>18</sup>O and N<sub>2</sub>H<sup>+</sup> 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<sup>18</sup>O and N<sub>2</sub>H<sup>+</sup>, our models do not succeed in reproducing the N<sub>2</sub>D<sup>+</sup> 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<sub>2</sub>D<sup>+</sup> is similar to that of N<sub>2</sub>H<sup>+</sup> and that the deuteration fraction [N<sub>2</sub>D<sup>+</sup>]/[N<sub>2</sub>H<sup>+</sup>] ~ 0.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<sub>2</sub>D<sup>+</sup>]/[N<sub>2</sub>H<sup>+</sup>] 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<sub>2</sub>H<sup>+</sup> and change the [N<sub>2</sub>D<sup>+</sup>]/[N<sub>2</sub>H<sup>+</sup>] ratio. In fact, the large linewidths observed (~1.5 km s<sup>-1</sup>) 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<sub>3</sub><sup>+</sup> isotopologues, and lowers the D-fractionation, as discussed by Caselli et al. (2002). (ii) The ortho-to-para H<sub>2</sub> 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<sub>2</sub> leads to a lower D-fractionation, given that the more energetic ortho-H<sub>2</sub> can more easily drive the proton-deuteron exchange reactions (e.g. H<sub>3</sub><sup>+</sup> + HD → H<sub>2</sub>D<sup>+</sup> + H<sub>2</sub>) backward and reduce the H<sub>2</sub>D<sup>+</sup>/H<sub>3</sub><sup>+</sup> abundance ratio (Gerlich et al. 2002). (iii) A higher ionization rate, $\zeta$, which may be due to the presence of X-rays. A larger $\zeta$ implies a larger electron fraction and thus a higher dissociative recombination rates for molecular ions, including H<sub>3</sub><sup>+</sup> 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<sub>2</sub>H<sup>+</sup> deuteration in a sample of representative IM Class 0 protostars. Our results can be summarized as follows.

- We observed the millimeter lines of C<sup>18</sup>O, C<sup>17</sup>O, N<sub>2</sub>H<sup>+</sup>, and N<sub>2</sub>D<sup>+</sup> 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<sub>2</sub>H<sup>+</sup> 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<sub>2</sub>, is low, below a few 0.001 in all Class I sources. There is, however, a wide dispersal in the values of R<sub>2</sub> 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<sub>2</sub>H<sup>+</sup> 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.
<figure><img src="image_16.png" /><figcaption>Fig. 15. The same as Fig. 11 for IC 1396 N and model 4.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. The same as Fig. 11 for CB 3 and model 3.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. The same as Fig. 11 for NGC 7129–FIRS 2 and model 3.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. The same as Fig. 11 for Serpens–Firs 1 and model 3.</figcaption></figure>
- 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, Tₑv. Our observational data are better fit with values of Tₑv ~ 20–25 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 $L_⊙$ to 84 $L_⊙$ (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, $n = 5 \times 10^5$ cm⁻³, $T_d = 13$ K) and the values assumed at $r = 15 000$ AU, which are $n = 1 \times 10^4$ cm⁻³ and $T_d = 10$ 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_0$, that is expected to be close to the canonical value ($X_0 = 2 \pm 2 \times 10^{-7}$) for radii lower than a given radius, $R_0$ and (ii) an abundance, $X_1$, that is expected to be <$X_0$ for larger radii. Thus

defined $R_0$, is the evaporation radius of CO. The values of $X_0$, $R_0$, and $X_1$ 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_1$ and $X_2$, 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 ~1.2 × 10⁻⁷; (ii) an intermediate layer with a high value for the C¹⁸O depletion, $f_D$ ~ 20, 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⁻¹⁰ 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_2D^+) = 1.9 \times 10^{-11}$ in the molecular cloud, and obtained an upper limit of $X(N_2D^+) < 2.0 \times 10^{-12}$ 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 $M_⊙$ and 80–120 $L_⊙$, 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 6 × 10⁻⁸ 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 $f_D$, ~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 ~7 × 10⁻¹² 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 6 × 10⁻⁸ 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 $f_D$ is ~5, and it would occur close to R = 5500 AU. The N₂H⁺ emission was better fit with a two-step function, $X_{N_2H^+} = 4.2 \times 10^{-10}$ for R < 10 000 AU, $1.8 \times 10^{-10}$ for R > 15 000 AU, and $1.0 \times 10^{-11}$ 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 <7 × 10⁻¹². The lower value of $f_D$ 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⁻³) 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 1.3 × 10⁻⁷ for radii less than 25 000 AU, decreases to <3×10⁻⁸ ($f_D$ > 3), and then increases to 9.0×10⁻⁷ 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 $X_{N_2H^+} = 6.5 \times 10^{-10}$. The N₂D⁺ emission is fit with a constant abundance of $X_{N_2D^+} = 3.0 \times 10^{-11}$, 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 $f_D = 3$, 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 ~4.5 × 10⁻¹¹, and the deuterium fractionation, ~0.1, is the typical value found in prestellar cores.

### References

* Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197

* Bechis, K. P., Harvey, P. M., Campbell, M. F., & Hoffmann, W. F. 1978, ApJ, 226, 439

* Beltrán, M. T., Girart, J. M., Estalella, R., Ho, P. T. P., & Palau, A. 2002, ApJ, 573, 246

* Beltrán, M. T., Girart, J. M., Estalella, R., & Ho, P. T. P. 2004, A&A, 426, 941

* Bergin, E. A., Maret, S., van der Tak, F. F. S., et al. 2006, ApJ, 645, 369

* Bisschop, S. E., Fraser, H. J., Öberg, K. I., van Dishoeck, E. F., & Schlemmer, S. 2006, A&A, 449, 1297

* Boogert, A. C. A., Pontoppidan, K. M., Knez, C., et al. 2008, ApJ, 678, 985

* Casali, M. M., Eiroa, C., & Duncan, W. D. 1993, A&A, 275, 195

* Caselli, P., Hartquist, T. W., & Havnes, O. 1997, A&A, 322, 296
Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002, ApJ, 565, 344

Caselli, P., Vastel, C., Ceccarelli, C., et al. 2008, A&A, 492, 703

Caux, E., Ceccarelli, C., Castets, A., et al. 1999, A&A, 347, L1

Ceccarelli, C., & Dominik, C. 2005, A&A, 440, 583

Cesaroni, R., Felli, M., Jenness, T., et al. 1999, A&A, 345, 949

Chini, R., Ward-Thompson, D., Kirk, J. M., et al. 2001, A&A, 369, 155

Codella, C., & Bachiller, R. 1999, A&A, 350, 659

Codella, C., Bachiller, R., Nisini, B., Saraceno, P., & Testi, L. 2001, A&A, 376, 271

Collings, M. P., Dever, J. W., Fraser, H. J., McCoustra, M. R. S., & Williams, D. A. 2003, ApJ, 583, 1058

Crapsi, A., Caselli, P., Walmsley, C. M., et al. 2005, ApJ, 619, 379

Crimier, N., Ceccarelli, C., Lefloch, B., & Faure, A. 2009, A&A, 506, 1229

Crimier, N., et al. 2010, A&A, accepted

Daniel, F., Cernicharo, J., Roueff, E., Gerin, M., & Dubernet, M. L. 2007, ApJ, 667, 980

de Gregorio-Monsalvo, I., Gómez, J. F., Suárez, O., et al. 2006, AJ, 132, 2584

Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803

Edwards, S., & Snell, R. L. 1983, ApJ, 270, 605

Eiroa, C., Palacios, J., & Casali, M. M. 1998, A&A, 335, 243

Emprechtinger, M., Caselli, P., Volgenau, N. H., Stutzki, J., & Wiedner, M. C. 2009, A&A, 493, 89

Felli, M., Palagi, F., & Tofani, G. 1992, A&A, 255, 293

Flower, D. R., Pineau des Forêts, G., & Walmsley, C. M. 2004, A&A, 427, 887

Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2005, A&A, 436, 933

Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2006, A&A, 449, 621

Froebrich, D., Smith, M. D., Hodapp, K.-W., & Eislöffel, J. 2003, MNRAS, 346, 163

Fuente, A., Martin-Pintado, J., Bachiller, R., Neri, R., & Palla, F. 1998, A&A, 334, 253

Fuente, A., Neri, R., Martín-Pintado, J., et al. 2001, A&A, 366, 873

Fuente, A., Martin-Pintado, J., Bachiller, R., Rodriguez-Franco, A., & Palla, F. 2002, A&A, 387, 977

Fuente, A., Neri, R., & Caselli, P. 2005a, A&A, 444, 481

Fuente, A., Rizzo, J. R., Caselli, P., Bachiller, R., & Henkel, C. 2005b, A&A, 433, 535

Fuente, A., Ceccarelli, C., Neri, R., et al. 2007, A&A, 468, L37

Fuente, A., Castro-Carrizo, A., Alonso-Albi, T., et al. 2009, A&A, 507, 1475

Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590

Fukui, Y., Iwata, T., Mizuno, A., Ogawa, H., & Takaba, H. 1989, Nature, 342, 161

Gerlich, D., Herbst, E., & Roueff, E. 2002, Planet. Space Sci., 50, 1275

Guesten, R., & Marcaide, J. M. 1986, A&A, 164, 342

Guillet, V., Jones, A. P., & Pineau Des Forêts, G. 2009, A&A, 497, 145

Hartigan, P., & Lada, C. J. 1985, ApJS, 59, 383

Harvey, P. M., Wilking, B. A., & Joy, M. 1984, ApJ, 278, 156

Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167

Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 263, 589

Henning, T., Burkert, A., Launhardt, R., Leinert, C., & Stecklum, B. 1998, A&A, 336, 565

Hily-Blant, P., Walmsley, M., Pineau Des forêts, G., & Flower, D. 2010, A&A, 513, A41

Huard, T. L., Weintraub, D. A., & Sandell, G. 2000, A&A, 362, 635

Hugo, E., Asvany, O., & Schlemmer, S. 2009, J. Chem. Phys., 130, 164302

Hurt, R. L., & Barsony, M. 1996, ApJ, 460, L45

Jones, A. P., Tielens, A. G. G. M., & Hollenbach, D. J. 1996, ApJ, 469, 740

Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2005, A&A, 437, 501

Larsson, B., Liseau, R., Men'shchikov, A. B., et al. 2000, A&A, 363, 253

Launhardt, R., & Henning, T. 1997, A&A, 326, 329

Lee, H.-H., Bettens, R. P. A., & Herbst, E. 1996, A&ASS, 119, 111

Lefloch, B., Eisloffel, J., & Lazareff, B. 1996, A&A, 313, L17

Maret, S., Ceccarelli, C., Caux, E., et al. 2004, A&A, 416, 577

Maret, S., Bergin, E. A., & Lada, C. J. 2006, Nature, 442, 425

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

Mckee, C. F. 1989, ApJ, 345, 782

Melnick, G. J., & Bergin, E. A. 2005, Adv. Space Res., 36, 1027

Minchin, N. R., White, G. J., & Ward-Thompson, D. 1995, A&A, 301, 894

Miranda, L. F., Eiroa, C., & Gomez de Castro, A. I. 1993, A&A, 271, 564

Neri, R., Fuente, A., Ceccarelli, C., et al. 2007, A&A, 468, L33

Noriega-Crespo, A., Moro-Martin, A., Carey, S., et al. 2005, ESA Spec. Publ., 577, 453

Öberg, K. I., van Broekhuizen, F., Fraser, H. J., et al. 2005, ApJ, 621, L33

Oliveira, C. M., Hébrard, G., Howk, J. C., et al. 2003, ApJ, 587, 235

Pagani, L., Vastel, C., Hugo, E., et al. 2009, A&A, 494, 623

Pagani, L., Salez, M., & Wannier, P. G. 1992, A&A, 258, 479

Patel, N. A., Greenhill, L. J., Herrnstein, J., et al. 2000, ApJ, 538, 268

Rodríguez, L. F., Moran, J. M., Ho, P. T. P., & Gottlieb, E. W. 1980, ApJ, 235, 845

Rodríguez, L. F., Curiel, S., Moran, J. M., et al. 1989, ApJ, 346, L85

Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369

Serabyn, E., Guesten, R., & Mundy, L. 1993, ApJ, 404, 247

Sipilä, O., Hugo, E., Harju, J., et al. 2010, A&A, 509, A98

Sugitani, K., Matsuo, H., Nakano, M., Tamura, M., & Ogura, K. 2000, AJ, 119, 323

Tofani, G., Felli, M., Taylor, G. B., & Hunter, T. R. 1995, A&AS, 112, 299

Umebayashi, T., & Nakano, T. 1990, MNRAS, 243, 103

Van der Tak, F. F. S., & van Dishoeck, E. F. 2000, A&A, 358, L79

Van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

Vastel, C., Caselli, P., Ceccarelli, C., et al. 2006, ApJ, 645, 1198

Wang, Y., Evans, N. J., II, Zhou, S., & Clemens, D. P. 1995, ApJ, 454, 217

Wilking, B., Mundy, L., McMullin, J., Hezel, T., & Keene, J. 1993, AJ, 106, 250

Wilson, T. L., & Matteucci, F. 1992, A&AR, 4, 1

Wilson, T. L., & Rood, R. 1994, A&ARv, 32, 191

Yun, J. L., & Clemens, D. P. 1994, AJ, 108, 612