The chemical structure of the Class 0 protostellar envelope NGC 1333 IRAS 4A

TL;DR

This study investigates the chemical structure of the Class 0 protostellar envelope NGC 1333 IRAS 4A, focusing on outflow influence, age estimation, and comparison with high-mass protostellar chemistry using multi-line observations and modeling.

Contribution

It provides detailed empirical molecular abundance profiles and compares them with chemical models, highlighting the roles of UV radiation and temperature in chemical differentiation.

Findings

Outflow activity influences deuteration and molecular abundances.

Estimated age of IRAS 4A is greater than 4×10^4 years.

Chemical abundances are lower than in high-mass envelopes, likely due to temperature differences.

Abstract

It is not well known what drives the chemistry of a protostellar envelope, in particular the role of the stellar mass and the outflows on its chemical enrichment. We study the chemical structure of NGC 1333 IRAS 4A in order to (i) investigate the influence of the outflows on the chemistry, (ii) constrain the age of our object, (iii) compare it with a typical high-mass protostellar envelope. In our analysis we use JCMT line mapping and HIFI pointed spectra. To study the influence of the outflow on the degree of deuteration, we compare JCMT maps of HCO+ and DCO+ with non-LTE (RADEX) models in a region that spatially covers the outflow activity of IRAS 4A. To study the envelope chemistry, we derive empirical molecular abundance profiles for the observed species using the radiative transfer code (RATRAN) and adopting a 1D dust density/temperature profile from the literature. We compare our best-fit observed abundance profiles with the predictions from the time dependent gas grain chemical code (ALCHEMIC). The CO, HCN, HNC and CN abundance require an enhanced UV field which points towards an outflow cavity. The abundances (wrt H2) are 1 to 2 orders of magnitude lower than those observed in the high mass protostellar envelope (AFGL 2591), while they are found to be similar within factors of a few with respect to CO. Differences in UV radiation may be responsible for such chemical differentiation, but temperature differences seem a more plausible explanation. The CH3OH modeled abundance profile points towards an age of > 4x10^4 yrs for IRAS 4A. The spatial distribution of H2D+ differs from that of other deuterated species, indicating an origin from a foreground colder layer (<20 K). The observed abundances can be explained by passive heating towards the high mass protostellar envelope, while the presence of UV cavity channels become more important toward the low mass protostellar envelope.

Tables

Table 1. Table 1: References for the collisional rate coefficients of the species modeled in this paper.

Species	Authors
CO	Yang et al. (2010)
HCO+	Flower (1999)
N2H+	Flower (1999)
HCN	Dumouchel et al. (2010)
CS	Lique et al. (2006); Lique & Kłos (2011)
HNC	Dumouchel et al. (2010)
H2CO	Wiesenfeld & Faure (2013)
CH3OH	Rabli & Flower (2010)
C2H	Müller et al. (2005)
CN	Lique et al. (2006); Lique & Kłos (2011)
OCS	Green & Chapman (1978)
SO2	Green (1995)
H2CS	Wiesenfeld & Faure (2013)

Table 2. Table 2: Constant empirical abundances estimated with RATRAN for the envelope of NGC 1333 IRAS 4A based on HIFI and JCMT observations. The table contains also the abundances of AFGL 2591 (outer envelope) for direct comparison.

Molecule	Abundance wrt H2		Abundance wrt CO
	IRAS 4A	AFGL 2591	IRAS 4A	AFGL 2591
		(outer)
CO	3 $\times$ ${10}^{-5}$	2 $\times$ ${10}^{-4}$	1	1
HCO+	1 $\times$ ${10}^{-9}$	3 $\times$ ${10}^{-8}$	3.3 $\times$ ${10}^{-5}$	1.5 $\times$ ${10}^{-4}$
N2H+	8 $\times$ ${10}^{-11}$	8 $\times$ ${10}^{-10}$	2.7 $\times$ ${10}^{-6}$	4 $\times$ ${10}^{-6}$
HCN	3 $\times$ ${10}^{-10}$	5 $\times$ ${10}^{-7}$	1 $\times$ ${10}^{-4}$	2.5 $\times$ ${10}^{-3}$
CS	3 $\times$ ${10}^{-9}$	4 $\times$ ${10}^{-8}$	1 $\times$ ${10}^{-4}$	2 $\times$ ${10}^{-4}$
HNC	8 $\times$ ${10}^{-11}$	1 $\times$ ${10}^{-8}$	2.7 $\times$ ${10}^{-6}$	5 $\times$ ${10}^{-5}$
H2CO	4 $\times$ ${10}^{-10}$	1 $\times$ ${10}^{-8}$	1.3 $\times$ ${10}^{-5}$	5 $\times$ ${10}^{-5}$
CH3OH	1 $\times$ ${10}^{-8}$	8 $\times$ ${10}^{-8}$	3.3 $\times$ ${10}^{-4}$	4 $\times$ ${10}^{-4}$
DCO+	1 $\times$ ${10}^{-11}$	…	3.3 $\times$ ${10}^{-7}$	…
C2H	6 $\times$ ${10}^{-10}$	8 $\times$ ${10}^{-8}$	3 $\times$ ${10}^{-5}$	4 $\times$ ${10}^{-4}$
CN	1 $\times$ ${10}^{-10}$	1 $\times$ ${10}^{-9}$	3.3 $\times$ ${10}^{-6}$	5 $\times$ ${10}^{-6}$
OCS	6 $\times$ ${10}^{-9}$	4 $\times$ ${10}^{-8}$	2 $\times$ ${10}^{-4}$	2 $\times$ ${10}^{-4}$
SO2	4 $\times$ ${10}^{-10}$	5 $\times$ ${10}^{-9}$	1.3 $\times$ ${10}^{-5}$	2.5 $\times$ ${10}^{-5}$
HDCO	5 $\times$ ${10}^{-11}$	…	1.7 $\times$ ${10}^{-6}$	…
H2CS	4 $\times$ ${10}^{-10}$	4 $\times$ ${10}^{-9}$	1.3 $\times$ ${10}^{-5}$	2 $\times$ ${10}^{-5}$
D2CO	5 $\times$ ${10}^{-11}$	…	1.7 $\times$ ${10}^{-6}$	…
H2S	4 $\times$ ${10}^{-10}$	4 $\times$ ${10}^{-9}$	1.3 $\times$ ${10}^{-5}$	2 $\times$ ${10}^{-5}$

Table 3. Table 3: As Table  2 , but for jump and drop abundance profiles.

	IRAS 4A		AFGL 2591		Jump/
Molecule	XIN	XOUT	XIN	XOUT	Drop at
CO	1 $\times$ ${10}^{-4}$	3 $\times$ ${10}^{-6}$	2 $\times$ ${10}^{-4}$ (const.)		25 K
H2CO	4 $\times$ ${10}^{-8}$	4 $\times$ ${10}^{-10}$	1 $\times$ ${10}^{-11}$	1 $\times$ ${10}^{-8}$	100 K
CH3OH	1 $\times$ ${10}^{-8}$	7 $\times$ ${10}^{-10}$	8 $\times$ ${10}^{-7}$	8 $\times$ ${10}^{-8}$	100 K

Table 4. Table 4: Initial abundances for chemical modeling.

Species	Abundance
“Low metals” (LM), $C/O=0.44$
o-H2	$3.75\times$ 10-1
p-H2	$1.25\times$ 10-1
HD	$1.55\times$ 10-5
He	$9.75\times$ 10-2
C	$7.86\times$ 10-5
N	$2.47\times$ 10-5
O	$1.80\times$ 10-4
S	$9.14\times$ 10-8
Si	$9.74\times$ 10-9
Na	$2.25\times$ 10-9
Mg	$1.09\times$ 10-8
Fe	$2.74\times$ 10-9
P	$2.16\times$ 10-10
Cl	$1.00\times$ 10-9

Table 5. Table 5: The top 25 of the initially abundant molecules for the NGC 1333 IRAS4A chemical modeling.

Species	Abundances
PSC-LM model with $C/O=0.44$
p-H2	$3.77\times$ 10-1
o-H2	$1.23\times$ 10-1
He	$9.75\times$ 10-2
H	$5.25\times$ 10-4
H2O∗	$5.53\times$ 10-5
CO∗	$4.05\times$ 10-5
CO	$3.26\times$ 10-5
O2	$1.79\times$ 10-5
HD	$1.52\times$ 10-5
N2	$7.39\times$ 10-6
NH${}_{}{}^{\ast }{}_3{}^{}$	$5.64\times$ 10-6
O	$5.59\times$ 10-6
O${}_{}{}^{\ast }{}_2{}^{}$	$4.12\times$ 10-6
CH${}_{}{}^{\ast }{}_4{}^{}$	$3.64\times$ 10-6
N${}_{}{}^{\ast }{}_2{}^{}$	$1.76\times$ 10-6
H∗	$6.03\times$ 10-7
C3H${}_{}{}^{\ast }{}_2{}^{}$	$4.48\times$ 10-7
OH	$3.43\times$ 10-7
H2O	$2.79\times$ 10-7
HNO∗	$2.40\times$ 10-7
NO	$2.22\times$ 10-7
N	$1.36\times$ 10-7
HDO∗	$1.35\times$ 10-7
CO2	$1.32\times$ 10-7
CO${}_{}{}^{\ast }{}_2{}^{}$	$1.19\times$ 10-7

Table 6. Table 6: Molecular Line Detections for NGC 1333-IRAS 4A - HIFI data

					Component 1			Component 2
Molecule	Transition	Rest Frequency	Eup	Vlsr	FWHM	Tmb	Vlsr	FWHM	Tmb
		(MHz)	(K)	(km s-1)	(km s-1)	(K)	(km s-1)	(km s-1)	(K)
CO	6-5	691473.08	116.2	$7.5\pm 0.1$	$8.04\pm 0.07$	$2.78\pm 0.01$	$6.6\pm 0.6$	$0.8\pm 0.2$	$5.2\pm 1.2$
13CO	6-5	661067.28	111.0	$6.7\pm 0.1$	$10.7\pm 0.4$	$0.19\pm 0.01$	$6.93\pm 0.01$	$1.60\pm 0.04$	$0.74\pm 0.01$
13CO	7-6	771184.12	148.1	$7.4\pm 0.3$	$9.0\pm 1.0$	$0.17\pm 0.02$	$6.78\pm 0.03$	$1.4\pm 0.1$	$0.48\pm 0.03$
C18O	6-5	658553.28	110.6	$6.8\pm 0.1$	$2.4\pm 0.3$	$0.14\pm 0.02$	$6.81\pm 0.05$	$0.6\pm 0.2$	$0.11\pm 0.03$
HCO+	8-7	713341.23	154.1	$6.1\pm 0.6$	$7.0\pm 2.0$	$0.17\pm 0.04$	$7.07\pm 0.04$	$1.9\pm 0.1$	$0.68\pm 0.05$
N2H+	7-6	652095.57	125.2	$7.7\pm 0.2$	$5.3\pm 0.3$	$0.05\pm 0.08$	$7.5\pm 0.1$	$1.5\pm 0.4$	$0.15\pm 0.08$
N2H+	8-7	745209.87	161.0	$7.9\pm 0.1$	$2.3\pm 0.3$	$0.13\pm 0.02$	$7.73\pm 0.05$	$0.5\pm 0.10$	$0.11\pm 0.02$
H2O	2(1,1)-2(0,2)	752033.14	136.9	$8.6\pm 0.6$	$30.2\pm 0.7$	$0.34\pm 0.01$	$0.2\pm 0.2$	$11.4\pm 0.9$	$0.24\pm 0.02$
H2S	2(1,2)-1(0,1)	736034.10	55.1	$6.9\pm 0.2$	$4.1\pm 0.2$	$0.13\pm 0.01$	$6.6\pm 0.1$	$0.5\pm 0.1$	$0.4\pm 0.1$
					Single component
HCN	8-7	708877.21	153.1	$6.4\pm 0.5$	$8.9\pm 0.5$	$0.10\pm 0.01$
HCN	9-8	797433.66	191.4	$4.7\pm 0.6$	$6.4\pm 0.6$	$0.09\pm 0.01$
CS	13-12	636531.84	213.9	$6.8\pm 0.5$	$6.2\pm 0.5$	$0.05\pm 0.02$
CS	14-13	685434.76	246.8	$4.5\pm 0.5$	$6.2\pm 0.5$	$0.05\pm 0.2$
CS	15-14	734324.00	282	$5.5\pm 0.6$	$6.1\pm 0.6$	$0.04\pm 0.01$
H2CO	9(1,9)-8(1,8)	631702.81	163.6	$7.0\pm 0.2$	$2.7\pm 0.3$	$0.18\pm 0.02$
H2CO	9(0,9)-8(0,8)	647081.73	156.2	$6.5\pm 0.2$	$3.8\pm 0.5$	$0.09\pm 0.01$
H2CO	9(5,5)-8(5,4)	655212.10	451.2	$6.2\pm 0.2$	$2.8\pm 0.6$	$0.08\pm 0.02$
H2CO	9(1,8)-8(1,7)	674809.78	174.0	$6.3\pm 0.2$	$4.6\pm 0.6$	$0.13\pm 0.01$
H2CO	10(1,10)-9(1,9)	701370.46	197.3	$6.9\pm 0.1$	$4.7\pm 0.2$	$0.19\pm 0.01$
CH3OH	5(-2,4)-4(-1,4)	665442.45	60.7	$6.5\pm 0.8$	$5.2\pm 0.8$	$0.05\pm 0.01$
CH3OH	4(2,3)-3(1,2)–	673745.93	60.9	$6.9\pm 0.4$	$4.2\pm 0.4$	$0.07\pm 0.03$
CH3OH	8(1,8)-7(0,7)++	674990.42	97.4	$5.9\pm 0.3$	$5.0\pm 0.3$	$0.09\pm 0.01$
CH3OH	4(2,2)-3(1,3)++	678785.45	60.9	$6.6\pm 0.5$	$6.1\pm 0.5$	$0.05\pm 0.01$
CH3OH	9(1,9)-8(0,8)++	719664.88	118.1	$6.9\pm 0.2$	$3.7\pm 0.3$	$0.09\pm 0.01$
CH3OH	5(2,3)-4(1,4)++	728862.52	72.5	$5.2\pm 0.5$	$5.0\pm 0.5$	$0.07\pm 0.01$
CH3OH	5(3,2)-4(2,2)	772453.80	82.5	$6.4\pm 0.6$	$5.6\pm 0.6$	$0.05\pm 0.01$
CH3OH	6(2,4)-5(1,5)++	779380.51	86.5	$6.9\pm 0.2$	$3.8\pm 0.2$	$0.15\pm 0.01$
					Absorption component
CO	6-5	691473.08	116.2	$7.5\pm 0.3$	$0.7\pm 0.1$	$-4.7\pm 2.1$
H2S	2(1,2)-1(0,1)	736034.10	55.1	$7.6\pm 0.3$	$0.9\pm 0.1$	$-0.3\pm 0.1$

Equations

$$\frac{N_u}{g_u}=x\frac{\int T_{mb}dV}{v{\mu }^{2}S}=\frac{N_T}{Q(T_{rot})}{e}^{-\frac{E_u}{T_{rot}}},$$

$$ln\frac{N_u}{g_u}=ln\frac{N_{T,thin}}{Q(T_{rot})}-\frac{E_u}{k{T}_{ex}}-ln(C_{\tau })+ln(f),$$

Full text

*Aims. *We study the chemical structure of NGC 1333 IRAS 4A in order to (i) investigate the influence of the outflows on the chemistry, (ii) constrain the age of our object, (iii) compare it with a typical high-mass protostellar envelope.

Methods. In our analysis we use JCMT line mapping and HIFI pointed spectra. To study the influence of the outflow on the degree of deuteration, we compare JCMT maps of HCO+ and DCO*+* with non-LTE (RADEX) models in a region that spatially covers the outflow activity of IRAS 4A. To study the envelope chemistry, we derive empirical molecular abundance profiles for the observed species using the radiative transfer code (RATRAN) and adopting a 1D dust density/temperature profile from the literature. We compare our best-fit observed abundance profiles with the predictions from the time dependent gas grain chemical code (ALCHEMIC).

Results. The CO, HCN, HNC and CN abundance require an enhanced UV field which points towards an outflow cavity. The abundances (wrt H2) are 1 to 2 orders of magnitude lower than those observed in the high mass protostellar envelope (AFGL 2591), while they are found to be similar within factors of a few with respect to CO. Differences in UV radiation may be responsible for such chemical differentiation, but temperature differences seem a more plausible explanation. The CH3OH modeled abundance profile points towards an age of $>$ 4$\times$104 yrs for IRAS 4A. The spatial distribution of H2D*+* differs from that of other deuterated species, indicating an origin from a foreground colder layer ($<$ 20 K).

*Conclusions. *The observed abundances can be explained by passive heating towards the high mass protostellar envelope, while the presence of UV cavity channels become more important toward the low-mass protostellar envelope.

1 Introduction

During low-mass ($<$ 2 M*⊙*) star formation a rotating cloud of gas and dust collapses under gravitational forces. The central protostar increases in mass through the accretion disk that surrounds it. The main mechanisms that retard the gravitational collapse are the thermal pressure, magnetic fields, and turbulence (Hennebelle & Motte 2009; Evans 2011; Luhman 2012; Tan 2015). Turbulence can be enriched by energetic outflows from young stellar objects (YSOs) which may further trigger star formation in nearby gas (Quillen et al. 2005).

Molecular outflows are prominent during the earliest stages of star formation, especially when collimated jets are driven in the youngest (103–104 years) embedded protostars (Arce et al. 2007). Class 0 protostars are still in their main accretion phase and they also drive the most powerful outflows. The impact of the ejected material on the surrounding cloud causes shock fronts. These lead to changes in the chemical composition and the enhancement of the abundance of several species in the surroundings. Fontani et al. (2014) have found enhancement of HDCO/H2CO towards the shock location of a Class 0 object, L1157 mm (d$=$250 pc), reporting a deuterated molecule as a shock tracer for the first time.

The strong outflow activity and winds that YSOs produce result in high-velocity gas, but also in the evacuation of regions near the protostar. Such cavities have been previously seen as a “hole” in the 1.3 mm continuum emission (e.g., near NGC 1333 IRAS4 and SVS13; Lefloch et al. 1998). UV radiation from the protostellar system (mainly due to accretion) is expected to play a crucial role in such environments (Stäuber et al. 2004; Visser et al. 2012).

During the cold and dense pre-collapse phase, molecular complexity increases by rapid ion-molecule gas-phase reactions followed by gradual freeze-out and build up of ices (H2O, CO, NH3) and surface processes. While collapsing, the radiation that comes from the forming protostar heats the inner parts of the envelope making surface radicals mobile and highly reactive. Later, these freshly formed complex ices thermally desorb, further boosting rich chemical processes in the gas and creating a “hot corino” (e.g., Ceccarelli 2004; Garrod & Herbst 2006; Ceccarelli 2008). The hot corinos refer to inner regions ($<$200 au) with an increase of the temperature above 100 K, as a result of passive heating from the protostar. An analog characterizes the high-mass protostars. Previous independent studies on high- and low-mass protostellar envelopes suggest that there is an abundance difference of a few orders of magnitude in several species (e.g., H2CO, CH3OH van der Tak et al. 2000; Maret et al. 2004, 2005). In addition, Herbst & van Dishoeck (2009) present complex organic molecule abundances relative to CH3OH that are similar (within factors of a few) for low- and high-mass YSOs.

In this work we are interested in a direct abundance comparison among high-mass and low-mass protostellar envelopes with the use of datasets from the same instruments and similar methodology. In particular, we aim to answer a) what the chemical structure of low mass protostellar envelopes is and how it compares to high mass protostellar envelopes, b) how does the temperature profile affect the abundance profile of several species in the inner envelope (“hot corino”), the freeze–out zones, and the outer parts of protostellar envelopes, c) how do the outflows influence the chemistry of the surroundings of a protostellar envelope and what is the role of outflow cavities in the observed abundances and finally, d) if the deuteration can be used as an outflow/shock tracer. For this purpose we use the low-mass protostellar envelope IRAS 4A, which appears as the brightest (sub)–mm continuum object in NGC 1333 IRAS 4 region and is classified as a Class 0 object (André et al. 1993). IRAS 4A is a prototypical well studied Class 0 object and of great interest as it is among the first (Mathieu 1994) proto–binary systems ever detected. NGC 1333 is one of the nearest (D=235 pc; Hirota et al. 2008) and youngest ($<$ 1 Myr; Gutermuth et al. 2008) star forming regions. IRAS 4A is a binary system, consisting of two deeply embedded Class 0 YSOs with a separation of 1.$\arcsec$8 (420 au at a distance of 235 pc). The binary nature of IRAS 4A was first observed in 0.84 mm CSO–JCMT interferometric high–resolution submillimeter continuum observations (Lay et al. 1995) and resolved at millimeter wavelengths using the BIMA array by Looney et al. (2000). They were also found to share a common circumbinary envelope (Looney et al. 2003).

In addition, a spectral line and continuum survey using the SMA interferometer was performed by Jørgensen et al. (2007) where inverse P–Cygni 13CO 2–1 line profiles have been found. These profiles indicate infall motions, which are also characteristic of the Class 0 stage. Di Francesco et al. (2001) reported inverse P–Cygni profiles in CS and H2CO, tracing high-density gas as observed with IRAM Plateau de Bure.

IRAS 4A has been suggested to have a “hot corino” (Maret et al. 2004). Multitransition observations of species such as H2CO and CH3OH towards 4A revealed abundance enhancements in the warmest inner regions ($>$100 K) by up to 2 orders of magnitude (Maret et al. 2004, 2005). The same abundance enhancement can also occur in outflows on larger scales as a result of ice mantle sputtering in shocks (Bachiller & Pérez Gutiérrez 1997; Tafalla et al. 2000). Mantle sputtering is thought to play a role when outflow speeds reach about 10 km s*-1* and is independent of gas density. In faster shocks with speeds as high as 20–25 km s*-1* the mantles vaporize completely (Guillet et al. 2011) while grain core sputtering is still inefficient. H2CO and CH3OH have been found to trace outflow activity of IRAS 4A (Jørgensen et al. 2007; Koumpia et al. 2016), which makes it difficult to distinguish between a “hot corino” chemistry and the enhancement due to shocks caused by protostellar outflows.

The highly collimated outflows from IRAS 4A have been mapped in several CO transitions (Knee & Sandell 2000; Jørgensen et al. 2007; Yıldız et al. 2012). IRAS 4A shows two bipolar outflows, one with a N–S orientation and the other with P.A. $\sim$45*∘. The usual interpretation of the CO and SiO outflow has been that a N–S component arises in A2* which becomes bent to a P.A. $\sim$45*∘* angle at short distance to the North. There is no evidence for a flow from A1, the brighter of the millimeter sources. Marvel et al. (2008) have used water masers to trace the small-scale motions of the IRAS 4A outflows, with somewhat puzzling results. One interpretation has been that there may be a third component of the system quite close to A2, which drives the outflow.

Our study mainly aims to a) constrain the chemical structure of the protostellar envelope of IRAS 4A and compare it with a high-mass protostellar envelope, b) investigate the presence of deuteration towards the outflow, and c) investigate the influence of the temperature over outflow activity to the observed abundance profiles. In this article we present HIFI and JCMT observations of a range of chemically diverse species towards IRAS 4A. To study the outflow chemistry, we estimate the excitation temperature and column density of H2CO in the envelope and outflow using population diagrams. We proceed with modeling our observed maps with RADEX in order to determine the DCO*+/HCO+* abundance ratio and test the enhancement observed by Fontani et al. (2014). To study the envelope chemistry, we apply the 1D physical model determined by Kristensen et al. (2012) that takes into account the temperature and density gradients of the envelope and we run RATRAN in order to fit our observations. We apply constant and jump-like abundance profiles. In addition, we run chemical models that similarly take into account the physical structure of the protostellar envelope and apply the abundance profiles we determine to our RATRAN models.

We compare this low–mass case with the high–mass case of AFGL 2591 for which a similar analysis was performed by Kaźmierczak-Barthel et al. (2015). Lastly, we try to constrain the chemical age of the IRAS 4A from our time-dependent chemical models.

2 Observations and data reduction

2.1 HIFI: observations and reduction

Observations of the YSO NGC 1333 IRAS 4A (RA=03h29m10.3s, Dec=$+$31o13′31″[J2000]) were made with the HIFI instrument at the Herschel Space Observatory, as part of the CHESS guaranteed time key programme (Chemical Herschel Survey of Star Forming Regions; Ceccarelli et al. 2010). Full spectral scans were made in bands 2a and 2b, covering the spectral ranges 626.01–721.48 GHz (415.81–479.23 $\mu$m) and 714.02–800.90 GHz (374.58–420.16 $\mu$m), respectively, at a resolution of 1.1 MHz (0.41–0.53 km s*-1*). In this frequency range, the Herschel/HIFI Half Power Beam Width (HPBW) is 26.5–33.9″(Roelfsema et al. 2012). The spectral scan was performed in Dual Beam Switch (DBS) mode with a normal chop frequency of 0.17 Hz. Each sky frequency was covered four times to facilitate the double sideband deconvolution process. The instrument stability settings were calculated in HSPOT by fixing the minimum and maximum goal resolutions to 1.1 MHz and setting the 1 GHz reference option without continuum optimization.

The observations were processed with the pipeline at the Herschel Science Center with HIPE 7.1.0 and retrieved from the Herschel Science Archive. Further post–pipeline level 2 processing was done in HIPE 8.0. Since we are interested in full spectral coverage, only the WBS spectra were considered, although many HRS spectra in narrower ranges were also available. Spectral regions affected by ‘spurs’ and not automatically detected by the HIFI pipeline were flagged out and ignored. Polynomial baselines were subtracted using the FitBaseline task by masking the lines interactively. The overlapping sidebands were deconvolved with the doDeconvolution task in HIPE by applying the default settings. The observed line intensities in units of antenna temperature were corrected for loss in the sidelobes by converting them to main beam temperatures using the main beam efficiency of $\eta_{mb}$ (Beff/Feff) of 0.75 (Roelfsema et al. 2012).

2.2 James Clerk Maxwell Telescope (JCMT): observations and reduction

These observations are a part of the JCMT Spectral Legacy Survey (SLS; Plume et al. 2007). The Auto–Correlation Spectral Imaging System (ACSIS) was used at the James Clerk Maxwell Telescope (JCMT111The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the Netherlands Organization for Scientific Research, and the National Research Council of Canada.) on Mauna Kea, Hawaii.

We have 2′$\times$2′ maps from the HARP–B instrument, which provides high velocity resolution (1 MHz, $\sim$ 1 km s*-1*). These observations cover the frequency window between 330 and 373 GHz. The angular resolution of JCMT is \sim 15$$\arcsec at 345 GHz, which is equivalent to $\sim 3500$ au at the distance of NGC 1333 IRAS 4 (Choi et al. 2004). The beam efficiency is 0.63 (Buckle et al. 2009).

Details regarding the reduction and line detections of this dataset can be found in Koumpia et al. (2016).

3 Observational results

3.1 Line detections

The Single Side Band (SSB) H and V–polarization HIFI spectra were searched independently for line detections. We consider as safe detections the signals detected in both polarizations with $>$ 3$\times$RMS (RMS $\sim$ 0.01 K–0.04 K) after averaging, and with widths of at least two channels $>$ 0.9 km s*-1* (single channel; 1.1 MHz, $\sim$0.47 km s*-1*). Table 6 presents the line list with secure detections. The detected lines were identified by producing single temperature LTE models in CASSIS222CASSIS has been developed by IRAP-UPS/CNRS (http://cassis.irap.omp.eu) of the species detected in the JCMT Spectral Legacy Survey by Koumpia et al. (2016). The species that were identified following this process are: CO, 13CO, C18O, CS, HCN, HCO*+, N2H+, H2CO, CH3OH, H2S and H2O. We inspected our HIFI spectra for lines of SO, SO2*, SiO, HNC, and H2O isotopologs, but found nothing. This is not surprising considering that the predicted intensity of the transitions in the HIFI range, assuming a beam filling of unity, is comparable with the measured RMS for HNC and a few orders of magnitude lower than the measured RMS for SO, SO2 and SiO (RMS$<$0.03 K). The transitions of these latter species are expected to show such weak lines in the 626–801 GHz regime mainly because of their a) high critical densities ($>$108 cm*-3*) and/or b) low Einstein coefficients ($<$10*-5* s*-1*) compared to the observed transition in the JCMT regime (330–373 GHz).

The JCMT observations have been described in more detail in Koumpia et al. (2016), where the detected species towards IRAS 4A are also presented. In addition to the species detected in the HIFI survey and excluding H2S and H2O, we detect the deuterated species DCO*+, HDCO, D2CO, the S–bearing species SO, SO2*, OCS and finally SiO, HNC, CN, C2H. Finally, we detect o–H2D*+* emission in the vicinity of IRAS 4A, but not directly towards the source. The overall rms noise level in the JCMT data ranges between 0.005 and 0.05 K at a velocity resolution of 0.9 km s*-1*.

3.2 Line profiles

Examples of the line profiles from the species detected with HIFI and JCMT are plotted in Figures 1–3. Large line profile variations are observed. The CO 6–5 line shows narrow emission (FWHM$\sim$0.8 km s*-1*) and absorption peaks accompanied by prominent wings extending over 25 km s*-1* towards blue- and red-shifted velocities. The absorption is absent in the other CO isotopologs, and the wings are much weaker. On the other hand, the H2O 211–202 line shows only very broad emission, very similar to the CO 6–5 wings (Figure 1).

Visually, the detected lines can be divided into three groups. These groups are characterized by line profiles that are dominated by a) absorption features such as CO 6–5 (Figure 1) and H2S 212–101, b) broad emission (FWHM$=$4–11 km s*-1*; Figure 2) such as 13CO 6–5, CS 13–12, H2O 211–202, HCN 8–7, CH3OH 532–422 and H2CO 918–817, and c) narrow emission (FWHM$<$4 km s*-1*; Figure 3) such as C18O 6–5, H2S 212–101, HCO*+* 8–7 and N2H*+* 8–7.

To quantify the difference between the line profiles, they were decomposed with a simultaneous fit of three Gaussians: broad emission, narrow emission, and narrow absorption components (Table 6). Our observations show that there is considerable variation in the widths of the broad emission line components. With median FWHM values of 4.6 and 5.1 km s*-1*, the H2CO and CH3OH lines are narrower than the H2O line (11 km s*-1*). The H2O line itself is in fact asymmetric, and can be fitted with an additional, narrower Gaussian (VLSR $=$ 0.55 km s*-1*, FWHM$=$3.86 km s*-1*). Finally, it is worth mentioning that the peak positions of the narrow N2H*+* lines are similar to the absorption component seen in CO 6–5 and H2S ($\sim$7.5 km s*-1*). The existence of such variation in the shape of the lines can be a result of the different regions that these lines trace. A broad component is indicative of outflow activity, a narrow component arises from dynamically quiescent gas (i.e., envelope) and the absorption is due to infall motions or the presence of foreground material. This foreground material can be a foreground cloud not eventually bended to the source itself, but it can also be the external and colder part of the envelope which absorbs the emission of the inner and warmer parts of this envelope.

4 Chemistry of the IRAS 4A outflow

4.1 LTE column densities

The protostellar envelope of NGC 1333 IRAS 4A has previously been studied by various authors (e.g., Maret et al. 2004, 2005), which is not the case for its associated outflow. In this section we aim to estimate the column densities and excitation temperatures from the observed emission in the envelope and outflow of IRAS 4A and compare them. In our approach we also add the observed HIFI transitions that were not available in previous studies.

A widely used method to estimate the column density of a molecule is the population diagram (Goldsmith & Langer 1999). When LTE applies, the Tex equals the gas kinetic temperature, otherwise it provides only a lower limit.

The column density of the upper state Nu and the rotational temperature Trot ($=$Tex) are determined by:

[TABLE]

where x $=$ 8.591$\times$1037 8$\pi$k/hc3, Nu the column density of the upper energy level (cm*-2*), gu the degeneracy of the upper energy level, Tmb the main beam temperature (K), dV the velocity range (km s*-1*), $\nu$ the frequency (Hz), $\mu$ the dipole moment, S the line strength, NT the total column density (cm*-2*), Trot the rotational temperature (K), Q the partition function and Eu the upper energy level. Plotting ln(Nu/gu) versus Eu/k results in a straight line with a slope of 1/Trot.

In non-LTE excitation the population of each level may be characterized by a different excitation temperature Tex ($\neq$Trot). The population diagrams can also take into account optical depth and beam effects due to different angular resolutions among the lines by using the modified equation:

[TABLE]

where Tex is the excitation temperature, Cτ$=$ $\tau$/(1-e*-τ*) is the optical depth correction factor, and f the beam dilution which is defined as the size of the telescope beam over the size of the emitting region and which is assumed equal for all lines. A more detailed description of the method and formulas used can be found in Goldsmith & Langer (1999).

The velocity information of our lines allows us to fit multiple Gaussians and separate, at first order, the outflow activity from the emission that comes from the dynamically quiescent envelope. For a more accurate determination of the temperature, species with many observed transitions are preferred. We perform this analysis for H2CO for which we have observed 12 transitions with JCMT (Koumpia et al. 2016) and HIFI (current study). In this approach we also include four extra transitions as observed using IRAM and presented in Maret et al. (2004).

Figures 4 and 5 present the population diagrams for the envelope and the outflow respectively, and the resulting excitation temperature ($\sim$rotational temperature), column density, optical depth and size of the emitting area. We find Tex\sim 35$$\pm0.5 K and N(H2CO) of 8.5{}^{+1.2}_{-0.8}$$\times 1014cm*-2* for the envelope and Tex\sim 41.5$$\pm4 K and N(H2CO) of 1.2{}^{+0.5}_{-0.2}$$\times 1015cm*-2* for the outflow. The error estimates have been computed during a $\chi^{2}$ minimization procedure. Maret et al. (2004) modeled only the envelope and found a factor of 2 lower column density and a factor of 1.2 lower Tex. We attribute these differences to the use of our HIFI transitions, which probe denser/warmer regions.

The excitation temperature provides only a lower limit for the gas kinetic temperature when the source is not in LTE. Our results point towards $\sim$20% lower Tex and $\sim$30% lower H2CO column density for the envelope compared to the outflow, which is significant given our error estimates ($\sim$10% and $\sim$14% respectively).

4.2 Deuterated species and outflows

When it comes to deuterated species in the field of star formation, the pre-stellar cores have been characterized as “deuterium fractionation factories” (Ceccarelli et al. 2014), due to the high deuteration degrees that have been observed towards them in various studies ($>$10% Bacmann et al. 2003; Crapsi et al. 2005b). Several studies have found a lower degree of deuteration towards Class 0 objects Crapsi et al. ($\sim$1%–10%; e.g. 2005a); Shah & Wootten ($\sim$1%–10%; e.g. 2001).

More recently, deuterated species have also been employed as shock tracers by Fontani et al. (2014). They found an enhancement of HDCO/H2CO ($\sim$10%) towards the eastern wall of the cavity excavated by the shock associated with the Class 0 object, L1157 mm. This is at least an order of magnitude larger than the HDCO/H2CO ratio of the surrounding material. We aim to study the distribution of the deuterated species in the region and examine the deuteration towards the outflow of NGC 1333 IRAS 4A, which is also a Class 0 object.

4.2.1 In search of ions and deuterated species towards the outflow

The outflow activity of IRAS 4A has been traced by many species. SiO is expected to be a tracer of outflow shocks due to sputtering of Si off dust grains (Schilke et al. 1997). Choi (2005) found it to trace the jet in the case of NGC 1333 IRAS 4A. They have presented a map of SiO 1–0 that traces the outflow activity of IRAS 4A and the spatial distribution of a narrow line component offset at $\sim$7.6 km s*-1*. Our dataset shows that SiO 8–7 has its primary peak at the shock position north of IRAS 4A (R1; Santangelo et al. 2014), while SO also emits significantly at that position (Figure 6). The same figure shows the integrated intensity map of C2H (core; from $+$5 to $+$9 km s*-1*) which was found to trace the envelopes of the three protostars. Its three peaks are in alignment with the continuum peaks as observed by Sandell & Knee (2001).

Figures 8 and 9 present the spatial distribution of the observed deuterated species D2CO, HDCO, DCO*+* and H2D*+* as overplotted with C2H. HDCO and D2CO trace mainly the protostellar envelopes and their distribution covers part of the outflow-shock area (white rectangle) north-northeast of IRAS 4A (R1; Santangelo et al. 2014). DCO*+* traces the envelopes but it also emits in a more extended area between the sources in the NW–SE direction. Interestingly, H2D*+* shows a very different spatial distribution compared to the other deuterated species. It does not follow the distribution of the envelopes and it mainly emits in the NW-SE direction in the space between the three sources. The spatial distribution of o-H2D*+* as presented in Figure 9 appears to be parallel with the narrow component of SiO as presented in Figure 1 by Choi (2005). Lastly, van der Tak et al. (2002) detected the rare ND3 species towards the DCO*+* peak (+23,-06 $\arcsec$ offset from IRAS 4A) which is about 15$\arcsec$ to the South of the H2D*+* peak (+38,+09 $\arcsec$ offset from IRAS 4A). Both ND3 and H2D*+* (Fig. 7) are characterized by narrow emission (1.0$<$FWHM$<$1.7 km s*-1*), indicative of quiescent gas or an outflow component perpendicular to the line of sight.

4.2.2 Origin of o-H2D*+* emission

H2D*+* is expected to arise from very cold gas. At very low temperatures ($<$20-30 K) the reaction H${}_{3}^{+}$ $+$ HD $\longleftrightarrow$ H2D*+* $+$ H2 $+$ $\Delta$E is not balanced by the backward process, increasing the abundance of H2D*+. In addition, the freeze out of CO and N2* normally boosts H${}_{3}^{+}$ (Roberts & Millar 2000), and increases the H2D*+* production rate (e.g., Bacmann et al. 2003; Caselli et al. 2003, 2008; Albertsson et al. 2013).

Given the “peculiar” spatial distribution of H2D*+, the logical follow-up is to investigate the spatial distribution of HCO+* and N2H*+* and their deuterated isotopologs. HCO*+* and N2H*+* are produced through the reactions: H${}_{3}^{+}$ $+$ CO $\Longrightarrow$ HCO*+* $+$ H2 and H${}_{3}^{+}$ $+$ N2 $\Longrightarrow$ N2H*+* $+$ H2. In dense and very cold environments both CO and N2 are expected to be depleted and thus so do HCO*+* and N2H*+. Given the fact that N2H+* is also destroyed by CO, the emission from the two species is usually spatially anti-correlated. Previous observations of N-bearing species (e.g. N2H*+) towards prestellar cores have shown depletion resistance compared to CO, and that N2* depletes at later times compared to CO (Bergin & Tafalla 2007; Pagani et al. 2012). This contradicts the expectation of a similar behavior of the two species in terms of freezing out and desorption mechanisms, which is based on the ratio of their binding energies ($\sim$1 under astrophysical conditions; Bisschop et al. 2006). More recent calculations and experiments on the topic are presented in Boogert et al. (2015).

Figure 10 presents the integrated intensity map of HCO*+* overplotted with H2D*+* and N2H*+* in contours. This figure shows that N2H*+* and HCO*+* have a similar spatial distribution, tracing mostly the protostellar envelopes and they emit significantly in almost half of the H2D*+* slab towards the N-NW axis.

The presence of N2H*+* and HCO*+* in part of the H2D*+* slab could be a result of the outflow-shock activity in the region that removes ices such as CO and N2 from grains. If we assume a single layer of gas, such activity should also make H2D*+* less abundant, but this is not what we observe. Previous studies including Choi et al. (2004) and Koumpia et al. (2016) provided evidence for a foreground cloud at $\sim$8 km s*-1* while IRAS 4A and IRAS 4B are part of a smaller embedded cloud at 6.7 km s*-1*. The H2D*+* line is shifted by $\sim$1.5 km s*-1* relative to the rest velocity of this cloud (Fig. 11). This could be an indication that the narrow component of SiO 1–0 ($\sim$7.6 km s*-1*) and the H2D*+* ($\sim$8 km s*-1*) originate from the same foreground layer at the offset velocity. Such narrow SiO emission has been discovered in more regions (e.g., G035.39–00.33) where it has been suggested that it originates from cold gas associated with a low-velocity shock (Duarte-Cabral et al. 2014). H2D*+* emission requires very cold conditions though and is unlikely to be associated with shock activity due to the presence of N2H*+* and HCO*+*. This emission probably originates from a colder layer in the foreground and the co-existence of the other two species is rather a projection effect.

4.2.3 Modeling the deuteration

Although high deuteration is usually connected with very cold environments, deuterated species have also been associated with shocks in a few studies (Lis et al. 2002; Fontani et al. 2014; Lis et al. 2016). We aim to model the distribution of the [DCO*+]/[HCO+*] ratio in the area covered by our JCMT maps, and examine the deuteration towards the outflow of NGC 1333 IRAS 4A.

In order to explore spatial variations in molecular D/H abundance ratios of the region, especially the [DCO*+]/[HCO+] ratio, we use the kinetic temperature map as derived by Koumpia et al. (2016) in the non-LTE radiative transfer program RADEX (van der Tak et al. 2007) after adopting a constant H2* density of 3$\times$105 cm*-3* as suggested in the same work. In addition to the optically thick HCO*+* (Koumpia et al. 2016) we use its optically thin isotopolog H13CO*+* which helps us derive an accurate [DCO*+]/[HCO+*].

We present the [DCO*+]/[HCO+] ratio in Figure 12. We find ratios varying between $\sim$0.01 and 0.07$\pm$0.016 around the protostellar envelopes and the surrounding gas covering the outflow towards the north of IRAS 4A. The error estimate reflects the average of the observational uncertainties, which takes its higher values as one moves to the edges of our maps where the signal from the lines becomes more comparable to the RMS. In particular, we find a $\sim$3 times higher [DCO+]/[HCO+] abundance ratio towards the N–NE part of the H2D+* slab compared to IRAS 4A and IRAS 4B while the N–NW part of the slab does not show an enrichment in deuteration. The shock position on the north though (R1) is characterized by values equal and up to two times higher in deuteration compared to IRAS 4A. The observed differences are significant given our error estimates. The observed [DCO*+]/[HCO+] abundance ratio is 3 orders of magnitude higher than the cosmic ratio [D]/[H] of 1.5$\times$10-5* (Linsky et al. 1995) and it reflects the strong deuteration occurring at those early embedded stages of star formation.

An enhancement of DCO*+* abundance requires an enhancement of H2D*+* abundance and the presence of CO in the gas phase, since it is produced through the reaction H2D*+* $+$ CO $\Longrightarrow$ DCO*+* $+$ H2. As we described in Section 4.2.2, H2D*+* is enhanced in very cold conditions($<$ 20–30 K). This is in agreement with the estimated kinetic temperature towards the shock position (20–30 K). Under the observed conditions CO freezes out though. To explain the observed DCO*+* abundance enhancement towards the shock position we conclude that this emission is dominated by deuterated species originally formed in the gas phase after the removal of CO from grains at low gas kinetic temperature ($<$ 20 K). Such release of CO into the gas phase at such low temperatures can occur during the passage of a shock wave.

5 Chemistry of the IRAS 4A envelope

5.1 RATRAN – model setup

In order to estimate molecular abundance profiles through the envelope of NGC 1333 IRAS 4A, we ran the Monte Carlo radiative transfer code RATRAN (Hogerheijde & van der Tak 2000) and produced synthetic line emission. RATRAN takes into account the physical structure of the source including temperature and density gradients, kinematics and continuum dust emission and absorption.

The line spectra of NGC 1333 IRAS 4A are generally broad ($>$5 km s*-1*) and the wings are very prominent in many species. The physical models are determined using continuum observations that characterize the protostellar envelope. Therefore we focus on the narrow component of the lines. For this purpose we perform a multi-Gaussian fit and use only the narrow component for modeling. Some of the lines show heavy self-absorption making their Gaussian fitting very inaccurate and thus, we chose not to include them in our models.

The H2O 211–202 line shows only a broad component (Figure 1) and we observe no isotopologs, making the identification of the envelope component unreliable. The situation for SiO is similar, thus we do not model these species. We preferably model the isotopologs when present, and otherwise model the narrow component of the main species when multiple Gaussian fitting is possible (e.g., H13CO*+, C17O and C18O, H13CN). We use fixed local ISM isotopic ratios of 12C/13C $=$ 60, 16O/17O $=$ 2000, and 16O/18O $=$ 560 (Wannier 1980; Wilson & Rood 1994; Wilson 1999). In addition we assume an ortho to para ratio of 3 for the two collisional partners o–H2* and p–H2. Dust continuum radiation is taken into account using the dust opacity OH5 taken from Ossenkopf & Henning (1994), which corresponds to dust grains with thin ice mantles.

We ran RATRAN applying the density and temperature radial profiles of NGC 1333 IRAS 4A (Figure 13) as determined by Kristensen et al. (2012). The density profile is assumed to be a power-law, n(r) $=$ n0 $\times$(r/r0)-p, with a best-fit index of p $=$ 1.8. The density and temperature profiles were derived from the best-fit dust model assuming that the gas is entirely molecular and using a mean molecular mass of 2.4 amu and a gas-to-dust ratio of 100. For our calculations we defined a grid of 19 spherical shells from rin $=$ 5$\times$1014 cm (33 au) up to rout $=$ 7.7$\times$1016 cm (5147 au) where the dust temperature is 250 K and 10 K and the H2 density is 3.05$\times$109cm*-3* and 3.5$\times$105cm*-3*, respectively. The original model extends further out as seen in Figure 13 but our adopted profiles better represent the size and outer temperature of the envelope of a protostar ($\sim$5000 au, 10 K). We also assumed thermal equilibrium between dust and gas at those high densities.

The papers from which the collisional rate coefficients of the main isotopologs with H2 were adopted are presented in Table 1. The collisional data of H2S are scaled from the files of ortho- and para-H2O as calculated by Dubernet et al. (2009); Daniel et al. (2010, 2011). For isotopologs and the deuterated species, the same collision data were used as for the main isotopolog.

We assumed a static envelope without infall or expansion and we used a range of abundances typically varying between 10*-7* and 10*-12* in order to constrain the abundance profile that best fits the observations. The turbulent line width was fixed to 1.9 km s*-1*, which is the average value we found for the narrow component of most species after fitting multiple Gaussians. Modeling the observed lines with a constant abundance is the only way when it comes to species that show only a few transitions, but might not always be a realistic approach. Several molecules have been suggested to be present in volatile ice mantles on dust grain surfaces at temperatures $<$ 20–110 K. The exact temperature depends on species and the surface composition (Bisschop et al. 2006; Herbst & van Dishoeck 2009). We chose to apply jump-like abundances at 100 K for H2CO and CH3OH for which more transitions are available and for which a constant abundance does not result in a good fit.

5.2 RATRAN – model results

The resulting abundances from the process described above are presented in Table 2 which compares this low-mass case with a high-mass case from the literature (AFGL 2591; Kaźmierczak-Barthel et al. 2015). We find observed abundance profiles for the low-mass protostellar envelope (NGC 1333 IRAS 4A) that are systematically 1-2 orders of magnitude lower than the high-mass protostellar envelope (AFGL 2591). Although, one can extract this information from previous studies on high and low-mass protostellar envelopes (e.g., H2CO, CH3OH van der Tak et al. 2000; Maret et al. 2004, 2005), this study is the first direct comparison among high-mass and low-mass protostellar envelopes with the use of datasets of the same instruments and similar methodology (radiative transfer and chemical models). We attribute the observed differences to the absence of a freeze-out zone (i.e., temperature differences) towards the high-mass protostellar envelope. Although AFGL 2591 is a much more distant object (3.3 kpc) than NGC 1333 IRAS 4A (235 pc), which could potentially affect our results, our models take the different distances into account.

Table 2 also presents the resulting abundances with respect to CO. In particular, we find a similarity in the observed abundances with respect to CO within a factor of a few among the low- and the high-mass protostellar envelope. Herbst & van Dishoeck (2009) present complex organic molecule abundances relative to CH3OH for low- and high-mass YSOs and find similar results. Our results suggest that gas phase CO/H2 measurements are essential for comparison among low- and high-mass protostellar envelopes in the future.

The observed data for CO, H2CO, and CH3OH do not fit well for constant abundances and are discussed below.

5.3 Drop models

We chose to work with CO isotopologs (C17O and C18O) because they are optically thin and have narrow line-widths and thus trace the quiescent envelope. The fact that we do not use the optically thick broad 13CO lines avoids systematic errors in our solution. Figure 14 shows the observed line profiles of C17O 3-2 and C18O 6-5 with the best fit abundance profile model overplotted after adopting a constant CO abundance of 3$\times$10*-5* (scaled accordingly for the isotopologs) and a drop model, where the CO abundance profile shows a drop in the freeze-out zone and increases again, for comparison. Using a constant abundance, it was impossible to reproduce the intensities of both C17O and C18O lines. In particular, the constant CO abundance reproduces the C18O 6–5 transition very well, but it overproduces the C17O 3-2 by a factor of $\sim$2. Having additional transitions would help to better constrain the abundance profile. The fact that the higher–J C18O 6–5 can be reproduced by a constant abundance profile indicates that this emission does not originate from the snowline of CO, which is characterized by a drop in the abundance.

Taking into account the freeze–out zone one would expect a drop of the CO abundance but this alone is also not able to reproduce the line intensities. Thus, we use a similar drop profile as suggested by Yıldız et al. (2010, 2012), in which the abundance of CO drops in the freeze-out zone but rises again in the outer envelope. For the evaporation temperature of CO we chose to use the lower limit from the laboratory, which is $\sim$25 K (Collings et al. 2003; Muñoz Caro et al. 2010; Luna et al. 2014). We find the best fit model is one with an inner abundance for T$>$25 K of 1$\times$10*-4* which drops to 3$\times$10*-6* at the coldest part of the envelope and rises again to the canonical value of 1$\times$10*-4* towards the outer envelope where external UV radiation becomes important. This is in good agreement with the additional C18O transitions from Yıldız et al. (2012) who found a CO abundance of 6$\times$10*-5* for the warmest part of the envelope (50$<$T$<$250 K) dropping to 3$\times$10*-6* for the coldest part ($<$30 K) where depletion is prominent and a jump up to 3$\times$10*-4* again for the outer envelope.

5.4 Jump models

The assumption of a constant abundance does not work well for all transitions of H2CO and CH3OH, indicating that a jump model is required. In jump models an abundance enhancement by up to a few orders of magnitude is applied in the inner warmest regions, where the temperature is sufficient to evaporate the grain mantles ($>$100 K).

Figure 15a shows the integrated intensities of the observed H2CO transitions using HIFI and JCMT towards IRAS 4A overplotted with the convolved ($\sim$15–33$\arcsec$) synthetic emission as calculated with RATRAN for a constant abundance of [H2CO] $=$ 4$\times$10*-10* and jump models at $\sim$100 K. Our jump model at $\sim$100 K results in XIN $=$ 4$\times$10*-8* and XOUT $=$ 4$\times$10*-10* and reproduces the lower-energy lines better but the higher-energy lines slightly worse compared to the constant abundance model. Maret et al. (2004) found H2CO abundance about a factor of two lower for the inner envelope (2$\times$10*-8*) and a factor of two lower for the outer envelope compared to our jump model. Differences up to 2 orders of magnitude between H2CO abundance of the inner and outer low-mass protostellar envelopes have been observed before (IRAS 16293–2422; Ceccarelli et al. 2000).

Figure 15b shows the integrated intensities of the observed CH3OH transitions using HIFI and JCMT towards IRAS 4A overplotted with the convolved ($\sim$15–33$\arcsec$) synthetic emission as calculated with RATRAN for [CH3OH] $=$ 1$\times$10*-8*. Maret et al. (2005) found an upper limit of 1$\times$10*-8* for the CH3OH abundance in the inner envelope and the same abundance as we do in the outer envelope.

The models we used are characterized by the same power law for the density profiles and the temperature profiles also follow a similar pattern. The main differences in the two methods are the additional HIFI transitions we use, the new collisional data, and the new physical model derived by Kristensen et al. (2012) including PACS observations, which were the smallest-scale used data. Maret et al. (2004, 2005) used the physical model derived by Jørgensen et al. (2002).

5.5 Deuteration

Using the abundances reported in Table 2 we find an HDCO/H2CO ratio of $\sim$10% which is comparable with the one found by Loinard et al. (2000, 2001) towards IRAS 16293–2422 ($\sim$10–15%) and the value reported for Orion KL ($\sim$15%) (Liu et al. 2011). We also find a very high D2CO over H2CO ratio ($\sim$10%) although this is based on only one transition. Very high values between $\sim$5–10% of the relative D2CO abundance have been found before towards the low-mass protostar IRAS 16293–2422 (Ceccarelli et al. 1998; Loinard et al. 2000), which is approximately 2 orders of magnitude higher than in Orion KL. The average gas kinetic temperature of IRAS 4A is $\sim$45 K (Koumpia et al. 2016) and thus such high deuteration cannot be explained if D2CO is only formed in the gas-phase. On the other hand, D2CO may have been enriched in the dust grains in low-mass protostellar environments during the cold, dense pre-collapse period followed by its evaporation thus it is probably formed via the grain surface reactions. Lastly, we find DCO*+/HCO+* $=$ 1% in agreement with what is reported by Loinard et al. (2000, 2001) for IRAS 16293–2422.

6 Chemical modeling

6.1 Pseudo time-dependent model

In this section we compare the empirical abundance profiles of the studied species obtained with RATRAN with those resulting from chemical models for the same adopted physical structure of NGC 1333 IRAS 4A. Our goal is to a) better understand the chemical processes that take place in NGC 1333 IRAS 4A and b) constrain its age.

To compare the observationally derived abundance profiles, we used the same 1D physical model as in Section 5.1 and the same gas-grain time-dependent chemical model as in Kaźmierczak-Barthel et al. (2015). This allows us to compare chemistry modeling results between a low-mass and a high-mass protostellar envelope. The chemical model is based on the chemical kinetics ALCHEMIC code of Semenov et al. (2010) and a gas-surface network with deuterium fractionation (Albertsson et al. 2013, 2014). The original, non-deuterated chemical network stems from the osu.2007 ratefile developed by the group of Eric Herbst 333http://web.archive.org/web/20081204232936/ (Garrod & Herbst 2006). The network is supplied with a set of approximately $1\,000$ high-temperature neutral-neutral reactions from Harada et al. (2010, 2012) and updated as of June 2013, using the KIDA database444http://kida.obs.u-bordeaux1.fr.

All H-bearing species reactions in this network were cloned by adding D, with the exception of molecules with the $-$OH functional group. Primal isotope exchange reactions for H${}_{3}^{+}$ as well as CH${}_{3}^{+}$ and C2H${}_{2}^{+}$ from Roberts & Millar (2000); Roberts et al. (2004); Roueff et al. (2005) were included. In cases where the position of the deuterium atom in a reactant or in a product was ambiguous, a statistical branching approach was used (for further details please consult Albertsson et al. 2013). In Albertsson et al. (2014) this deuterium network was further extended by adding ortho- and para- forms of H2, H${}_{2}^{+}$ and H${}_{3}^{+}$ isotopologs and the related nuclear spin-state exchange processes.

We adopted the cosmic ray ionization rate of $5\times 10^{-17}$ s*-1*, as derived by van der Tak & van Dishoeck (2000) and Indriolo et al. (2015). In addition, we ran chemical simulations with a higher value of $10^{-16}$ s*-1*. The photodissociation and photoionization rates in the model are adopted for a 1D slab model from van Dishoeck (2006). To mimic the presence of an outflow cavity in this source, we also considered several other models with enhanced UV irradiation. For that, we use a scaled interstellar UV radiation field in Draine (1978) units and moderate dust extinctions of 10, 3, 2, and 1 mag. The pre-computed photoionization and dissociation rates for a 1D plane-parallel slab model were used (see Eq. 5; Semenov et al. 2010). The self- and mutual-shielding of CO and H2 from external dissociating radiation was calculated as in Semenov & Wiebe (2011). Water self-shielding is neglected.

The grains are assumed to be uniform and spherical and made of amorphous olivine with a density of $3$ g cm*-3* and a radius of $0.1\,\mu$m. If grains were on average bigger than $0.1\,\mu$m it would slow down freeze-out which would increase the CO gas/ice ratio, especially for short timescales (1000 years; Figure 23). The same effect would be seen for all other species that can freeze out. This would enhance the overall gas-phase molecular abundances with respect to our standard model.

Each grain provides $\approx 1.88\times 10^{6}$ surface sites (Biham et al. 2001) for freeze-out of gaseous molecules. We use the desorption energies Edes from (Garrod & Herbst 2006). We compute the diffusion energies of surface reactants by multiplying their binding (desorption) energies by 0.4. In the current literature, factors of 0.3–0.8 are commonly used, motivated by existing data on stable species (e.g., Cuppen et al. 2017), although a firm physical basis for such a scaling is lacking.

The gas-grain interactions include sticking of neutral species and electrons to dust grains with 100% probability and desorption of ices by thermal, cosmic ray, and UV-driven processes. In our model we use a single Edes for H2 $=$ 450 K, which properly describes H2 diffusion over the dust surface and the H2 binding to water or silicate surfaces. However, the H2–H2 binding energy is much lower (23 K). This means that as soon as there is a H2 monolayer on a surface, further freeze-out of H2 would be compensated by immediate desorption of H2 back to the gas phase. Our models cannot capture this process well and thus we simulate it by not allowing H2, HD, and D2 to stick to grains. There is still enough H2, HD, and D2 forming on grains though, and thus their surface abundances are not zero (See also; Hincelin et al. 2015; Wakelam et al. 2016).

The UV photodesorption yield of $10^{-3}$ was used (e.g., Öberg et al. 2009; Fayolle et al. 2011, 2013). This is consistent for CO and some other light species, but it drops to $\sim$10*-5* or lower for anything bigger (e.g., CH3OH; Martín-Doménech et al. 2016; Bertin et al. 2016). This is one of the limitations of our method. Photodissociation processes of solid species are taken from Garrod & Herbst (2006).

Surface recombination is assumed to proceed through the classical Langmuir-Hinshelwood mechanism (e.g., Hasegawa et al. 1992). We do not allow tunneling of surface species via the potential wells of the adjacent surface sites. To account for hydrogen tunneling through barriers of surface reactions, we have employed Eq. (6) from Hasegawa et al. (1992), which describes a tunneling probability through a rectangular barrier with a thickness of 1 Å.

For each surface recombination, we assume there is a 1% probability for the products to leave the grain due to the partial conversion of the reaction exothermicity into breaking the surface-adsorbate bond (Vasyunin & Herbst 2013). Following experimental studies on the formation of molecular hydrogen on amorphous dust grains by Katz et al. (1999), the standard rate equation approach to the surface chemistry is utilized. In addition, dissociative recombination and radiative neutralization of molecular ions on charged grains and grain re-charging are taken into account.

With this network and $10^{-5}$ relative and $10^{-20}$ absolute tolerances, the 1D IRAS 4A model takes about 1 minute using Core-i7 2.5 GHz CPU (OS X 10.11, gfortran 6-x64) to compute over $10^{6}$ years. This time span encompasses the likely age of this object.

6.2 Initial abundances

To set the initial abundances, we calculated the chemical evolution of a 0D molecular cloud with $n_{\rm H}=2\times 10^{4}$ cm*-3*, T $=10$ K, and $A_{\rm V}=10$ mag over 1 Myr (model “PSC-LM”). For that, the neutral “low metals” (LM) elemental abundances of Graedel et al. (1982); Agúndez & Wakelam (2013) were used, with the solar ${\rm C/O}=0.44$, initial ortho/para H2 of 3:1, hydrogen being fully in molecular form, and deuterium locked up in HD (see Table 5).

6.3 Error estimations

The problem of uncertainties of the calculated abundances is well known in chemical studies of various astrophysical environments, ranging from dark clouds to hot cores (see, e.g., Dobrijevic et al. 2003; Vasyunin et al. 2004; Wakelam et al. 2005; Vasyunin et al. 2008). The error budget of the theoretical abundances is determined by both the uncertainties in physical conditions in the object and, to a larger degree, by uncertainties in the adopted reaction rate coefficients and their barriers. Poorly known initial conditions for chemistry may also play a role here.

In order to estimate the chemical uncertainties rigorously, one needs to perform a Monte Carlo modeling by varying reaction rates within their error bars and re-calculating the chemical evolution of a given astrophysical environment. We do not attempt to perform such a detailed study and use the estimates from previous works.

In previous studies of chemical uncertainties it was found that the uncertainties are in general larger for bigger molecules as their evolution involves more reactions compared to simpler molecules. For simpler, key species such as CO and H2 involved in a limited cycle of reactions it is easier to derive the reaction rates with a high accuracy of $\sim 25\%$. In addition these species are formed in the gas, which is better known than the surface part. Consequently, their abundances are usually accurate within $10-30\%$ in modern astrochemical models. On the other hand, for other diatomic and triatomic species such as CN, HCO+, HCN, CCH, and so on, uncertainties are usually about a factor of 3-4 (see Vasyunin et al. 2004, 2008; Wakelam et al. 2010). The chemical uncertainties are higher for more complex molecules like methanol due to the fact that gas phase reactions are less known and the surface chemistry less well understood. These uncertainties can reach orders of a magnitude, with the factor of 10 being a likely lower limit.

Moreover, for S-bearing species, for which many reaction rates have not been properly measured or calculated or included in the networks, these intrinsic uncertainties and hence the uncertainties in their resulting abundances are higher, $\ga$ factor of 10 even for simple species such as SO, OCS, and SO2 (Loison et al. 2012). Also, the incompleteness of astrochemical networks with regard to the chemistry of Cl- and F-bearing molecules makes their calculated abundances rather unreliable.

In our study we assume that the uncertainties in the abundances of ortho- and para-H2 and CO are within a factor of 30%. For HCO+, H2CO, CN, N2H+, C2H, NO, OH, C, C+, O, CH, NH3, H2O, HCN, and HNC the uncertainties are within a factor of three, and for S-bearing species, CH3OH, and HCl they are within a factor of ten.

The result of this process is the abundance profiles from the species of interest. This means that we get the abundance of species over the radius adopted from the physical models for a range of timescales. In our study we do not use the abundance profiles of SO2, SO, and OCS from the chemical model because the chemical network for S-bearing species is too inaccurate with respect to these species and thus lacks predictive power.

6.4 Results

6.4.1 Standard approach

Figures 16–18 present the results of the standard chemical modeling compared to the observed abundances of NGC 1333 IRAS 4A from Sec. 5.1 and the high-mass protostellar envelope AFGL 2591 from Kaźmierczak-Barthel et al. (2015). The observed abundances for most species appear to be in agreement with the modeled abundances in the outer envelope while they are systematically 1 to 2 orders of magnitude lower than the high-mass protostellar envelope. In contrast, the predicted CN, HCN, and HNC abundances are 1 to 2 orders of magnitude higher than the observed values for the outer envelope. Our chemical models do not take into account shielding of CN by H2, as well as FUV scattering, which can be important. In addition our models use the CN photodissociation rates taken by van Dishoeck et al. (2006). A more recent study by el-Qadi & Stancil (2013) presents CN cross-sections with values several times smaller than those from van Dishoeck et al. (2006). From all the parameters, it seems that the strongest effect in those modeled abundances is due to the assumed FUV intensity and the C/O ratio. The C/O $=1.1$ gives X(CN) of 3$\times$10*-7* while X(HCN) $\sim$ X(HNC) $\sim$ 8$\times$10*-10*. If one assumes a strong FUV field of 104 of the ISRF UV with a modest extinction of Av$=$1 mag, the fit is much better: X(CN) $\sim$1–4$\times$10*-10*, X(HCN) \sim$$<$$10^{-11}, X(HNC) \sim$$> 2\times$$10^{-11}. Without any additional FUV one gets X(CN) $=$ 10*-14*, X(HCN) $\sim$ 2$\times$10*-9*, X(HNC) $\sim$ 7$\times$10*-11*. The standard model alone cannot explain the observed abundances for more than one of the modeled species, making the development of a more advanced model necessary.

6.4.2 The necessity of an outflow cavity

We observe a drop of only $\sim$2 orders of magnitude towards the snowline of CO, compared to $\sim$6 orders of magnitude predicted by our chemical models. A plausible explanation for such a discrepancy between our chemical model and the observations is the presence of the outflow, which is not accounted for in the 1D model (Bruderer et al. 2009). The way to approximately simulate the outflow, the UV-irradiated outflow walls, and the envelope in the 1D approximation is to add more UV radiation to the chemical model. For that, additional FUV components with intensities of 1 and 10 Draine’s units and moderate dust extinctions of 10 and 3 mag were considered. We find that only lower extinction influences the resulting CO abundance, increasing it by approximately 1 order of magnitude. The standard model without extra UV radiation and extinction of 10 mag produces the same abundance profiles.

Figures 19–21 present the results of the models with additional 10 Draine UV fields and for dust extinctions of 1, 2 and 3 mag. The extreme case of a 100 Draine UV field and A${}_{V}=1$ was also considered resulting in an almost constant CO abundance profile. Our observed CO abundance profile appears to be reproduced reasonably well by UV$=$ 10$\times$ISRF and A${}_{V}=$1 mag (green line; Figure 14). Although such a model makes the overall fit better for CO, HCO*+* and DCO*+, it does not really improve the fit for other species (e.g., CS, CN, C2H) and it actually makes the overall fit worse for H2CO, HDCO, D2CO and CH3*OH, by decreasing the abundance of the outer envelope by up to 4 orders of magnitude. Thus, we do not consider it to be our best–fit (or standard) model.

The modeled HCO+ abundances generally follow those of its parent molecule, CO, and show a strong decline for radii between $\sim$ 334 au and 5350 au, where the CO freeze-out zone is located. HCO+ abundances also drop strongly in the inner, dense and dark envelope region at $r\la$134 au, where the ionization degree drops due to fast recombination processes. Not surprisingly, the modeled N2H+ abundances also drop in the very inner envelope, like those of HCO+. In contrast to HCO+, N2H+ thrives in the CO freeze-out zone, where a key destruction reaction of N2H+ ions by CO molecules is no longer effective. In Figures 17 and 20 the DCO+ abundances are compared, and the modeled DCO+ profile follows the case of HCO+ and agrees with observations in two areas: 1) The inner envelope at $\sim$67–267 au and 2) the outer envelope with $r\ga$6685 au. In contrast to HCO+, DCO+ is slightly overproduced in the no-UV chemical model in the inner part of the envelope, but is well fitted by the model with additional UV due to the outflow cavity.

The poor fit and dependence on UV irradiation of CN, HCN, and HNC have already been discussed above. Their formation mainly proceeds via neutral-neutral gas-phase reactions involving light hydrocarbons like C2H and other N-bearing species (e.g., NO). Thus the no-UV model that fails to fit the CN, HCN, and HNC data is also not able to fit the C2H observed abundances.

The modeled CS abundance profile shows a poor fit to the data as well. As we also mentioned above, this is due to a general lack of predictive power of current astrochemical models for S-bearing species. Still, H2S modeled values are in good agreement with the observed data.

The observationally-driven H2CO abundances are only well reproduced in the outer envelope at $r\ga$2005 au, and are lower than the observed values by up to 3-5 orders of magnitude in the inner part. This is also likely caused by the same approximation of the outflow and UV-irradiated outflow cavity walls in our 1D chemical model as for CO. Alternatively, our observations may lack the necessary resolving power and sensitivity (that interferometers can provide) to uncover and unbiasedly constrain the underlying physical structure of the inner NGC 1333 IRAS 4A envelope, which comes out in the lack of agreement between the data and the chemical predictions. A similar behavior is shown by H2CO isotopologs and the chemically-related methanol molecule.

In Figures 19–21 we show the effect of including an additional UV component in the chemical model, as our attempt to represent the UV-irradiated outflow cavity material. As discussed above, additional UV radiation lowers the degree of CO depletion and brings us much closer to agreement between the observed and modeled CO abundances. The same effect is seen for our N-bearing species (CN, HCN, HNC). Unfortunately, for all other observed species (HCO+, N2H+, CS, H2CO isotopologs, CH3OH, C2H, S-bearing species) the modeled abundance profiles have poorer agreement with observations than in the standard model. The enhanced UV irradiation leads to overly rapid destruction of less tightly bound molecules than CO, and limits the efficiency of surface chemistry by desorbing ices too efficiently. The potential solution to such a chemical discrepancy is to perform chemical modeling using a more realistic 2D or 3D physical structure of the NGC 1333 IRAS 4A envelope, including the outflow and outflow cavity wall and performing UV radiative transfer.

6.4.3 Time dependence and different input parameters

Our models are time dependent so we also investigate the influence of different timescales on the produced chemical abundance profiles. Figure 22 in the appendix presents the time dependent abundance profiles of several species, demonstrating the insignificant influence of time in the short timescales that characterize a Class 0 object (104-105 yrs; e.g. Enoch et al. 2009). Our modeled methanol abundance though is in better agreement with the observed abundance for $\geq$ 4$\times$104 yrs (Figures 17, 22), while other species do not show significant abundance variation on these timescales and thus we cannot use them as additional constraints for age. We provide a lower limit to the age of IRAS 4A which is at least four times older than the one given by Maret et al. (2002) and potentially in agreement with the value of 9$\times$104 yrs given by Gonçalves et al. (2008) based on the morphology of the observed magnetic field. We should point out that the derived best–fit age of our object is dependent on the framework of model that we use. If in reality conditions are different (e.g., presence of a disk), or our chemical network misses some key reactions, the best-fit age value can vary significantly. A rigorous way to do such modeling would require running numerous models with varying temperature, density, CR ionization rate, and reaction rates, which would give a best-fit chemical age plus its error bars. A previous study towards young high-mass star-forming regions took these factors into account and found chemical ages that were characterized by uncertainties of a factor of 2–3 (Gerner et al. 2014).

To test the dependence of our results on the adopted physical conditions, we have run models with a twice higher cosmic ray ionization rate, a grain growth up to 0.5 $\mu$m, and different initial abundances for timescales between 103 and 106 yrs. To set the different initial abundances for the chemical modeling, we calculated the chemical evolution of a 0D model of an infrared dark cloud with nH $=$ 2$\times$ 105 cm*-3*, T $=$ 15 K, $\zeta$CRP $=$ 5$\times$10*-17* s*-1*, H2 OPR $=$ 3:1 and AV $=$ 10 mag over 1 Myr. The neutral “low metals” elemental abundances of Graedel et al. (1982) and Agúndez & Wakelam (2013) were used. The resulting abundance profiles of this process can be seen in Figure 23.

We find that the CO abundance profiles are not strongly affected by increasing the cosmic ray ionization rate. In contrast, the different initial abundances can cause a decrease of 1–2 orders of magnitude in the abundance profiles at the inner or outer envelope but do not significantly affect the CO abundance at the snow line. Only in combination with the short timescale of 103 yrs do we see 0.5–1 order of magnitude higher abundance at the snow line compared to the other timescales. This timescale is too short for a Class 0 object though. Lastly, the model with bigger grain sizes of 0.5 $\mu$m shows an increase in CO abundances by 1–2 orders of magnitude at the snow line. This is because the total dust surface area per unit gas volume is smaller in this model compared to the standard case of 0.1 $\mu$m grains, lowering the pace and hence the degree of the CO freeze-out. The effect is particularly dramatic at shorter chemical ages of $\sim 10^{3}$ years, where the difference in CO abundances between the 0.1 $\mu$m and 0.5 $\mu$m models is $\sim 3-4$ orders of magnitude. When grains are on average bigger than $0.1\,\mu$m it slows down freeze-out which would increase the CO gas/ice ratio especially for short timescales (1000 years; Figure 23). The same effect would be seen for all other species that can freeze out. This would enhance the overall gas-phase molecular abundances with respect to our standard model. Although IRAS 4A is still a very young object ($\sim$10000 yrs) its lifespan cannot be only 1000 years. Even if we were to assume such a short lifespan, it could not explain a grain growth up to $0.5\,\mu$m, but rather an ISM-like size up to 0.2–0.3 $\mu$m (Bianchi et al. 2003). Some studies of mm observations have shown that larger, mm-sized grains may exist in the envelopes of Class 0 protostars (Kwon et al. 2009; Chiang et al. 2012). However, Class 0 objects are highly embedded, making it difficult to eliminate the optically thick emission, which can cause an overestimation on the grain size. In addition our modeled abundances cannot explain the observed ones for timescales of 10000 yrs only by applying a grain growth.

In conclusion, the time-dependent models show no significant differences in abundance profiles for the timescales that are relevant to our object. Therefore, adopting an age of 10000 yrs does not introduce significant systematic errors. The same applies for the variation of the cosmic ray ionization rate. In contrast, we find that grain growth and different initial abundances have a more significant influence in the resulting abundance profiles, especially at short timescales (1000 yrs). Such short chemical age in combination with significant grain growth is rather unlikely though for a Class 0 object such as IRAS 4A.

7 Conclusions

We used HIFI and JCMT data to constrain the chemical structure of a low-mass protostellar envelope and compare it with a high-mass equivalent.

7.1 Results

• •

Constant or jump/drop-like empirical abundance profiles reproduce well our single-dish submillimeter observations. The abundance in the outer envelope is supported for most species by the predictions of our 1D time-dependent gas-grain chemical model.

• •

The presence of an outflow cavity with a strong UV field can explain the observed abundances of several species (e.g., CO, HCO*+, DCO+*) in the outer low-mass protostellar envelope, while passive heating is sufficient towards a high-mass protostellar envelope.

• •

The empirical abundance (with respect to H2) profiles for the low-mass protostellar envelope (NGC 1333 IRAS 4A) are systematically 1 to 2 orders of magnitude lower than the high-mass protostellar envelope (AFGL 2591). The overall warmer temperature profile of high-mass protostellar envelopes seems to drive this result. We find similar empirical abundances as soon as we estimate them with respect to CO.

• •

Population diagrams for H2CO indicate 20% lower Tex and 30% lower H2CO column density for the envelope compared to the outflow.

• •

High D2CO over H2CO ratio (10%) towards IRAS 4A points towards formation via grain surface reactions during the cold phase and not gas-phase chemistry. This is in agreement with what has been observed before towards IRAS 16293-2422.

• •

We find an enrichment of DCO*+* over HCO*+* ratio towards the shock position compared to the protostellar envelope. We attribute this result to the CO originally formed in the grains and later released into the gas phase (at T$<$20 K) during the passage of the shock wave.

• •

H2D*+* shows a different spatial distribution compared to the other deuterated species and a peak velocity at $\sim$8 km s*-1*. The most prominent explanation is that it is located in a different layer of gas than the clump that contains the protostars.

• •

The abundance profile of CH3OH provides a lower limit for the age of NGC 1333 IRAS 4A of 4$\times$104 yrs.

7.2 Discussion

The modeled chemical abundance profiles of the inner envelope are a few orders of magnitudes lower than the observed ones for all species. We find a decrease of about 2 orders of magnitude in the abundance of species with more observed transitions such as CO, H2CO and CH3OH in the CO depletion zone (outer envelope). Similar drops have been seen before by interferometric studies of other low-mass protostars (e.g., IRAS 16293–2422, L1448–C; Schöier et al. 2004; Jørgensen et al. 2005). Jørgensen et al. (2007) suggest that the emission of H2CO and CH3OH are related to the shocks caused by the protostellar outflows rather than being the result of compact low-mass “hot corinos”. More transitions of the other observed species will help us in the future to better constrain their observed profile.

Furthermore, we tried to simply simulate an outflow cavity by increasing the UV radiation that the observed species are exposed to. We found that this approach improved the fit among the theoretical abundance profiles and the observed for several species (e.g., CO, HCO*+*), thus a more detailed 2D/3D chemical modeling that takes into account disk structure and outflow cavities is expected to be more accurate. Such a model requires more transitions and interferometric (high sensitivity and high spatial resolution) observations that trace the inner region of the protostellar envelope.

Lastly, we attribute the observed abundance difference with respect to H2 among the low- and the high-mass protostellar envelope to the higher temperatures that characterize the high-mass case and the absence of a freeze-out zone. For safer comparison, further studies of the same nature between high-mass and low-mass protostellar envelopes are necessary. In particular, the similarity in the observed abundances with respect to CO, suggests that gas-phase CO/H2 measurements are essential, as all other species are off by within factors of a few as the one of CO.

Acknowledgements.

D.S. acknowledges support from the Heidelberg Institute of Theoretical Studies for the project Chemical kinetics models and visualization tools: Bridging biology and astronomy. The authors thank Lars Kristensen for useful discussions. We also thank Inga Kamp and Veronica Allen for careful reading who helped in improving the clarity of this paper.

Appendix A Line detections towards NGC 1333-IRAS 4A

Appendix B Chemical models for various input parameters