© ESO, 2012

1. Introduction

In the embedded phases of low-mass star formation (Stages 0 and I; Whitney et al. 2003; Robitaille et al. 2006), the material surrounding the protostar is exposed to energetic phenomena such as shocks and ultraviolet radiation (Shu et al. 1987; Spaans et al. 1995; Bachiller & Tafalla 1999; Arce et al. 2007). Observationally, this leads to much higher intensities in lines of rotationally excited carbon monoxide, water and other species than would be possible from gas that is only heated through collisions with warm dust in the bulk of the envelope. The ubiquity of surplus rotationally excited emission was first established through ground-based observations of the CO J = 6–5 transition (upper-level energy E_u/k = 120 K; Schuster et al. 1993; Hogerheijde et al. 1998). The Infrared Space Observatory (ISO) subsequently detected CO up to the 21–20 line (E_u/k = 1300 K), along with rotationally excited H₂O (up to 470 K) and OH (up to 620 K; Ceccarelli et al. 1999; Giannini et al. 1999, 2001; Nisini et al. 2002). The origin of this hot gas has been debated ever since. The shape of the CO 6–5 spectra – a narrow emission feature on top of a broader one – and the presence of a strong, narrow ¹³CO 6–5 emission feature point at a mixture of quiescent and shocked gas, with the quiescent gas likely heated by the protostellar UV field along the cavity walls (Spaans et al. 1995). Spatially resolved CO 6–5 maps, showing extended narrow emission, provide additional support for this scenario (van Kempen et al. 2009). However, it remains unclear at a quantitative level to what extent UV radiation and shocks each contribute to heating up the gas and powering the emission.

The Herschel Space Observatory (Pilbratt et al. 2010) offers a new set of tools to tackle this question. Amongst its first results are detections of highly rotationally excited CO (up to J = 38−37, E_u/k = 4100 K) and H₂O (up to 640 K) in the embedded low-mass protostars HH 46 and DK Cha (van Kempen et al. 2010a,b). A quantitative radiative transfer model of HH 46 was successful at explaining the CO line fluxes with a combination of UV-heated gas in the outflow cavity walls and small-scale C-type shocks along the walls (van Kempen et al. 2010b). However, the model was unable to disentangle the origin of the observed H₂O lines in a similar fashion. The goals of the current paper are (1) to present our model from van Kempen et al. in more detail; (2) to check how robust the model fit is for CO; (3) to derive a quantitative fit for H₂O; and (4) to extend the analysis to DK Cha and a third source (NGC 1333 IRAS2A) for which CO and H₂O data have since become available (Kristensen et al. 2010; Yıldız et al. 2010). The three test cases are all nearby protostars (d = 180–450 pc) and cover the evolutionary range from deeply embedded Stage 0 to late Stage I.

The question of the origin of the hot gas observed with Herschel is not limited to embedded low-mass young stellar objects (YSOs). Herschel observations of rotationally excited CO and/or H₂O include intermediate-mass YSOs (Fich et al. 2010; Johnstone et al. 2010), high-mass YSOs (Chavarría et al. 2010), circumstellar disks (Meeus et al. 2010; Sturm et al. 2010), the Orion Bar photon-dominated region (Habart et al. 2010), and even extragalactic sources like the ultraluminous infrared galaxy Mrk 231 (van der Werf et al. 2010). A quantitative explanation of the origin of the observed emission is still missing in most of these cases. Our model may help in constraining possible excitation mechanisms in environments other than low-mass YSOs. In particular, the treatment of how the gas is excited by UV fields and shocks can be adapted for different types of sources.

This work is part of the two complementary Herschel key programmes WISH and DIGIT. The primary goal of WISH, or “Water in star-forming regions with Herschel” (van Dishoeck et al. 2011), is to use H₂O as a probe of the physical and chemical conditions during star formation and to follow its abundance as a function of evolutionary phase. DIGIT, or “Dust, ice, and gas in time” (PI: N. Evans), aims to study the evolution of the gas during different stages of star formation, as well as that of the dust grains and their icy mantles. Observations of CO are an integral part of both programmes, because its simple chemistry and easy use in radiative transfer codes make CO an excellent tool to constrain physical conditions and processes. Indeed, the approach in the current paper is to fit the free parameters in the model to the CO data and then to assess how well the model reproduces the H₂O data. More generally, the model is ultimately intended as a tool in using the Herschel data as a probe of the energetics and dynamics of star formation, in particular how and where the young star deposits energy back into the parent molecular cloud. Future work in WISH and DIGIT will explore these questions in more detail.

A detailed description of the model appears in Sect. 2. The Herschel data and complementary ground-based data are summarised in Sect. 3. The model results for CO and H₂O are presented in Sect. 4 and discussed in Sect. 5. Conclusions are drawn in Sect. 6.

2. Model description

The far-infrared (far-IR) and sub-millimetre (sub-mm) molecular line emission from embedded low-mass protostars has traditionally been interpreted with purely spherical envelope models, but such models cannot reproduce the observed intensities of rotationally excited lines like CO 6–5 and higher (e.g., van Kempen et al. 2009). Our model goes a step further by introducing a bipolar outflow cavity with two additional heating mechanisms: UV irradiation of the gas in the cavity walls and small-scale C-type shocks along the walls (Fig. 1). Rotationally excited emission may also arise from the bipolar jet and the associated internal J-type shocks, from gas heated by X-rays, or from a circumstellar disk embedded within the envelope. However, our goal is not to create a single all-inclusive model of an embedded protostar. Rather, the main question in this work is what the Herschel observations can teach us about the passively heated envelope, the UV-heated cavity walls and the small-scale shocks. Investigating the quantitative contributions from other components will be left for future studies. Furthermore, the current focus is on the integrated line intensities only. With regards to the Herschel-PACS data, we are only interested in the emission from the central spaxel of the full 5 × 5 spaxel array. Future updates of the model will explore the spatial extent of the emission and could add dynamical processes to allow for the computation of line profiles.

2.1. Passively heated envelope

Low-J lines from CO and its isotopologues and similar low-excitation lines from other species are commonly modelled with a spherical envelope (e.g., Jørgensen et al. 2002; Schöier et al. 2002). A spherical envelope is also the starting point of our model, since the warm inner region can contribute to emission in the higher-J lines. The gas in the envelope is well coupled to the dust and is only heated passively or indirectly, i.e., the dust is heated by the protostellar luminosity and the gas attains the same temperature through collisions with the dust.

The density and temperature profiles of each source are constrained from continuum observations. Following the method of Jørgensen et al. (2002) and Schöier et al. (2002), the continuum emission for a grid of envelope models is calculated with the 1D continuum radiative transfer program DUSTY (Ivezić et al. 1999). The best-fit model is determined by comparing the emission to spectral energy distributions (SEDs) and sub-mm brightness profiles compiled by Froebrich (2005) and Di Francesco et al. (2008). The inner radius of the envelope is set at a dust temperature of 250 K (about 30 AU or 0.15′′ at 200 pc; Table 2). If any material is present on smaller scales, its contribution to the overall line emission would be too diluted in the  ~10′′Herschel beam to have a significant impact on the model results. The envelopes are optically thick at UV and visible wavelengths, so most of the stellar radiation is reprocessed by the dust. For a constant stellar luminosity, the assumed stellar temperature therefore does not have a significant effect on the resulting dust temperatures or the molecular line emission (Jørgensen et al. 2002). This was verified by running part of the DUSTY grid for stellar temperatures of both 5000 and 10 000 K.

The remaining three free parameters are fitted in a χ² procedure: the power-law index for the density profile (n ∝ r ^− p), the optical depth at 100 μm (τ₁₀₀), and the size of the envelope (Y ≡ r_out/r_in). For our line radiative transfer models, the envelopes are terminated at the point where the density reaches 10³ cm^-3 or the temperature reaches 10 K (Jørgensen et al. 2002). This point is denoted r_env. The best-fit values for the three test cases are listed in Table 2 in Sect. 3, and the associated density and temperature profiles are plotted in Fig. 5.

The Doppler broadening parameter b for the passive envelope models is fixed at 0.8 km s^-1, a typical value for the quiescent gas in low-mass YSOs as determined from optically thin C¹⁸O observations (Jørgensen et al. 2002). The ongoing collapse of the envelope onto the star and disk is simulated with a power-law infall velocity profile ν ∝ r ^− 1/2 (Shu 1977), scaled to 2.0 km s^-1 at the inner edge of the envelope. Details on how the abundances are obtained and how the molecular line intensities are calculated are provided in Sects. 2.2.2 and 2.2.3.

2.2. Photon-heated cavity walls

2.2.1. Density and temperature

The second model component is the UV-heated gas in the walls of the outflow cavities (red in Fig. 1; Spaans et al. 1995; van Kempen et al. 2009). For an ambient cavity density of a few 10³ cm^-3 (Neufeld et al. 2009), the total gas column from the central source to the point where the cavity intersects the outer edge of the envelope (a distance of  ~10¹⁷ cm) is a few 10²⁰ cm^-2. This corresponds to a UV optical depth of a few tenths, so UV photons produced near the protostar can freely reach the cavity walls. When they do, they produce a photon-dominated zone where the gas is heated to temperatures of a few hundred K.

The first step in modelling the emission from the UV-heated gas is to carve out a bipolar cavity from the spherical envelope, similar to the procedure followed by Bruderer et al. (2010). Mid-infrared Spitzer Space Telescope images of HH 46 and similar YSOs like L 1527 show ellipsoidal cavities, with a ratio between the semimajor and semiminor axes (a_cav and b_cav) of about 4 (Velusamy et al. 2007; Tobin et al. 2008). Large-scale CO 3–2 maps of NGC 1333 IRAS2A, tracing the swept-up outflow gas, show roughly the same shape (Sandell et al. 1994). DK Cha is oriented close to pole-on (Knee 1992), so the shape and size of its outflow are unknown.

Tests of our model show that the variation in the computed CO and H₂O line fluxes resulting from changes in the shape and size of the cavities is within the overall model uncertainties, consistent with the more detailed tests of Bruderer et al. (2010). Our approach is therefore to fix the ratio between the outflow axes (a_cav/b_cav) of all three sources to the value of 4.3 measured for HH 46 (Velusamy et al. 2007), and to scale the absolute length of the axes with the envelope radius. In other words, the ratios a_cav/r_env and b_cav/r_env of all three sources are fixed to the values of 2.1 and 0.50 derived for HH 46. The envelope radii and outflow axis lengths are summarised in Table 2.

The two cavities for each source are aligned along the z axis, with the ends meeting at z = 0 (Fig. 2). The ellipsoidal cavity wall in the upper right quadrant of a cylindrical coordinate system is defined as ${\mathit{R}}_{\mathrm{cav}}\mathrm{=}{\mathit{b}}_{\mathrm{cav}}\sqrt{\mathrm{1}\mathrm{-}{\left(\frac{{\mathit{z}}_{\mathrm{cav}}}{{\mathit{a}}_{\mathrm{cav}}}\mathrm{-}\mathrm{1}\right)}^{\mathrm{2}}}\mathit{.}$(1)The points along the wall can also be characterised by their distance r_cav to the protostar, ${\mathit{r}}_{\mathrm{cav}}\mathrm{=}\sqrt{{\mathit{R}}_{\mathrm{cav}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{cav}}^{\mathrm{2}}}\mathit{,}$(2)which is used for some figures in Sect. 4.

The only UV source included in our model is the protostar. UV photons can also be produced at the bow shock (Böhm et al. 1993; Raymond et al. 1997) or by small-scale shocks within the cavity or along the walls (Curiel et al. 1995; Saucedo et al. 2003; Walter et al. 2003), but these contributions are highly uncertain. The protostellar UV luminosity itself is also highly uncertain, and is therefore treated as a free parameter. Classical T Tauri stars have typical L_UV of 10^-3.5–10^-1.5 times the total stellar luminosity (Ingleby et al. 2011). Scaling up to the higher accretion rates found in embedded YSOs, anything from 0.01 to 1 L_⊙ would be a reasonable value for our model. The best-fit UV luminosities derived in Sect. 4.2 range from  < 0.03 to 0.08 L_⊙. An L_UV of 0.1 L_⊙ would give an unattenuated UV flux at 100 AU of G₀ ≈ 10⁴ relative to the mean interstellar radiation field of Habing (1968).

Due to the UV field impinging on the cavity walls, the gas is heated to higher temperatures than it attains in the passively heated, purely spherical envelope. The temperature can be calculated by solving the heating and cooling balance as function of the physical conditions. Many such models are available in the literature, but they show substantially different results, sometimes even at the qualitative level (Röllig et al. 2007). The consequences of these uncertainties for our model results are discussed extensively in Sect. 5.1. For the initial results in Sect. 4, the PDR surface temperature grid of Kaufman et al. (1999) is adopted as our standard prescription. This grid was calculated for densities up to 10⁷ cm^-3; for higher densities, it is extrapolated as 1/n(H₂) based on the expected dominant heating and cooling mechanisms (Sect. 5.1).

The Kaufman et al. (1999) grid as such can only be used at the cavity wall. At any depth into the envelope, the temperature needs to be corrected for attenuation of the UV field by dust. Röllig et al. (2007) plotted the temperature as function of the visual extinction, A_V, for several densities and unattenuated UV fluxes. In each case, the temperature decreases roughly as exp( − 0.6A_V). In order to use this depth dependence in our model, the amount of extinction towards any point (R,z) is calculated with the continuum radiative transfer code RADMC (Dullemond & Dominik 2004). The input stellar spectrum is a 10 000 K blackbody scaled to a given UV luminosity between 6 and 13.6 eV. This input spectrum is a reasonable approximation of the real spectrum emitted by a low-mass protostar with a UV excess due to accretion (e.g., van Zadelhoff et al. 2003). The local UV flux from RADMC, F_RADMC (also integrated from 6 to 13.6 eV), can be considered the product of the unattenuated flux and an extinction term, ${\mathit{F}}_{\mathrm{RADMC}}\mathrm{\left(}\mathit{R,z}\mathrm{\right)}\mathrm{=}{\mathit{F}}_{\mathrm{unatt}}\exp \mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathrm{UV}}\mathrm{)}\mathit{,}$(3)where F_unatt is the UV flux through a unit surface arising from a source at a distance $\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}}$. The UV optical depth, τ_UV, is effectively an average over the many possible paths a photon can follow from the source to the point (R,z). It is converted into a visual extinction via A_V = τ_UV/3.02 (Bohlin et al. 1978). Along with the UV field, RADMC also calculates the dust temperature along the cavity wall and throughout the envelope. This is used instead of the dust temperature from the 1D DUSTY models.

One caveat with RADMC is that it treats scattering of the radiation by dust as an isotropic process, while in reality grains are partially forward scattering (Li & Draine 2001). Tests with the radiative transfer code of van Zadelhoff et al. (2003), which does allow for anisotropic scattering, show flux differences of at most 20% between an isotropic scattering function and the scattering function of Li & Draine. This is within our overall model uncertainty.

The present goal is to reproduce the observed integrated line fluxes, not the detailed line profiles. Hence, the Doppler broadening and infall velocities for the UV-heated component are kept the same as in the passively heated component.

2.2.2. Abundances

Several thousand grid points are required to properly sample the envelope and the cavity walls in the 3D radiative transfer calculations (Sect. 2.2.3). It would be too time consuming to calculate the abundances of CO and H₂O for each individual point with a traditional chemical code. Instead, the abundances are obtained according to the method of Bruderer et al. (2009). A grid of chemical models is pre-computed for the appropriate range of densities, temperatures, unattenuated UV fluxes and extinctions. Abundances for each individual point in the envelope or the cavity wall are interpolated from this grid when constructing the source models for the radiative transfer code.

The basis of our chemical reaction network is the UMIST06 database (Woodall et al. 2007) as modified by Bruderer et al. (2009), except that X-ray chemistry is not included. The network allows all neutrals to freeze out onto the dust to account for the depletion of many gas-phase species in the colder parts of the envelope. Molecules can return to the gas phase through thermal desorption and photodesorption. Binding energies and photodesorption yields for CO and H₂O are taken from Fraser et al. (2001), Bisschop et al. (2006) and Öberg et al. (2007, 2009). Photodesorption occurs even in regions of high extinction because of a cosmic-ray-induced UV flux of about 10⁴ photons cm^-2 s^-1, equivalent to 10^-4 times the mean interstellar flux (Shen et al. 2004). In the gas phase, UV photons can dissociate or ionise many species; the relevant rates are calculated according to van Dishoeck et al. (2006) for a 10 000 K blackbody stellar spectrum, appropriate for a low-mass protostar with a UV excess. Self-shielding of H₂ and CO, which causes their dissociation rates to decrease more rapidly than those of other species, is included according to Draine & Bertoldi (1996) and Visser et al. (2009). The cosmic-ray ionisation rate of H₂ is set to 5 × 10^-17 s^-1 (Dalgarno 2006).

The grid of pre-computed abundances, from which abundances for the three sources are interpolated, is obtained by advancing the chemical reactions for a period of 10⁵ yr. As starting conditions, all carbon is locked up in CO, the remaining oxygen in H₂O, all nitrogen in N₂, and all remaining hydrogen in H₂. Reflecting the most likely formation pathways, CO and N₂ start in the gas phase, while H₂O starts frozen out on the dust. H₂ is predominantly formed on the grains, but because it rapidly evaporates even at 10 K, it starts in the gas phase. Elemental abundances are adopted from Aikawa et al. (2008) and Bruderer et al. (2009); the initial composition is summarised in Table 1.

2.2.3. Radiative transfer

Level populations and line intensities for the passively heated envelope and the UV-heated outflow cavity walls are calculated with the new full-3D radiative transfer code LIME (Brinch & Hogerheijde 2010). LIME uses an irregular grid with points sampled randomly but weighted according to a user-defined criterion such as density or temperature. This naturally results in a finely spaced grid in the regions most important for the astrophysical problem at hand.

The source geometry of our models introduces a substantial difficulty: the high-excitation CO and H₂O lines originate in a thin layer of hot gas along the cavity wall, so an accurate flux calculation requires a high grid point density inside and immediately outside this small region. Furthermore, the grid has to be able to resolve emission on scales ranging from  ~10 to  ~10⁴ AU. This cannot be accomplished with the default density or temperature weighting functions. Extensive convergence tests were carried out to find an alternative grid sampling scheme. The optimised grid consists of 30 000 points, picked randomly in log r between r_in and r_env. For every point with radial coordinate r, a point (R_cav,z_cav) is located on the cavity wall such that $\sqrt{{\mathit{R}}_{\mathrm{cav}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{cav}}^{\mathrm{2}}}\mathrm{=}\mathit{r}$. Let γ_cav be the angle between the point (R_cav,z_cav) and the midplane, i.e., tanγ_cav = z_cav/R_cav. The angular coordinate γ of each point is then determined as follows:

• 1.

  one third of the points are placed exactly at the cavity wall(γ = γ_cav);

• 2.

  one third of the points are picked with a random angle γ such that 0.9 < γ/γ_cav ≤ 1;

• 3.

  one third of the points are picked with a random angle γ such that 0 < γ/γ_cav ≤ 1 and then weighted by $\sqrt{\mathit{\gamma }\mathit{/}{\mathit{\gamma }}_{\mathrm{cav}}}$, which favours points near the midplane over points near the cavity wall.

The grid thus constructed is smoothed by five iterations of Lloyd’s algorithm (Lloyd 1982), which moves each point somewhat towards the centre of mass of its Voronoi cell. The purpose of this smoothing step is two-fold: (1) to turn some “needle-shaped” cells created by the random setup into more regular shapes (Brinch & Hogerheijde 2010); and (2) to push some grid points from the first sampling group into the cavity wall, so that the edge between the envelope and the cavity is well sampled on either side. The first two groups of points together cover the region of hot gas responsible for the more excited lines. The third group of points covers the bulk of the envelope, where the low-J emission originates.

With the grid established, LIME employs the standard two-step method for problems where local thermodynamic equilibrium (LTE) does not apply. The first step is an iterative procedure to calculate the rotational level populations based on the balance between radiative excitation, collisional excitation and collisional de-excitation. The second step is the ray tracing, i.e., calculating the emergent spectrum at a given distance and inclination. The collisional rate coefficients required for the first step are collected from the Leiden Atomic and Molecular Database¹ (LAMDA; Schöier et al. 2005), where the most recent data files for CO and H₂O are from Yang et al. (2010) and Faure et al. (2007), respectively. The ortho-to-para ratio for H₂ is taken to be thermalised, with a maximum value of 3. The ortho-to-para ratio of H₂O is fixed at 3. Dust is included in the model with a standard gas-to-dust ratio of 100 and the OH5 opacities from Ossenkopf & Henning (1994), appropriate for dust grains with a thin ice mantle.

Images are created at 0$\stackrel{\mathrm{″}}{.}$1 resolution (18–45 AU at the distances of the three sources; Table 2) to get a sufficiently fine sampling of the hot inner regions. In order to compare the model results to the observations, each image is convolved with a Gaussian beam profile. Beam sizes are calculated as the diffraction-limited beam for the wavelength and the telescope (Herschel, APEX or JCMT) at which the simulated line was observed. In the case of lines observed with Herschel-PACS, fluxes are extracted from the central $\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}\mathrm{\times }\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}$ to mimic the instrument’s central spaxel. For a given density, temperature and abundance structure, the estimated uncertainty on the CO line fluxes ranges from 30% for the 2–1 line to a factor of 2 for the high-J end of the ladder. For H₂O, it is a factor of 2 for all lines. The uncertainty is mainly due to the complex geometry and the small size of the emitting region for the high-J lines. The source models themselves lead to uncertainties of similar magnitude in the line fluxes, mostly because of the problems with the gas temperature discussed in Sect. 5.1.

Integrated intensities from APEX, the JCMT and Herschel-HIFI (${}^{\mathrm{\int }}\mathit{T}_{\mathrm{mb}}\mathrm{d}\nu$ in K km s^-1) and integrated fluxes from Herschel-PACS (${}^{\mathrm{\int }}\mathit{F}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }$ in W m^-2) can be mutually converted through $\mathrm{\int }{\mathit{F}}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }\mathrm{=}\frac{\mathrm{2}\mathit{k}\mathrm{\Omega }}{{\mathit{\lambda }}^{\mathrm{3}}}\mathrm{\int }{\mathit{T}}_{\mathrm{mb}}\mathrm{d}\nu \mathit{,}$(4)with Ω(λ) the solid angle subtended by the beam at the given wavelength. Denoting the full width at half maximum (FWHM) of the beam as Θ(λ), the solid angle for APEX, the JCMT and HIFI is (π/4ln2)Θ². The formula is more complicated for PACS, because the size of the spaxels has to be taken into account. For just the central spaxel, the solid angle is ${\mathrm{\Omega }}_{\mathrm{PACS}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}{\left[\frac{\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}}{\mathrm{erf}\mathrm{(}\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}\sqrt{\ln \mathrm{2}}\mathit{/}{\mathrm{\Theta }}^{\mathrm{)}}}\right]}^{\mathrm{2}}\mathit{,}$(5)with erf the error function and Θ(λ) the FWHM of the diffraction-limited Herschel beam.

2.3. Small-scale shocks

The third model component is a series of C-type shocks moving along the cavity walls. These shocks may be powered by the wind launched from close to the protostar or the surface of the inner disk, or by the expansion of the jet (e.g., Lada 1985; Shang et al. 2006; Arce & Sargent 2006; Santiago-García et al. 2009). Shocks are created when the wind hits the outflow cavity walls, heating the gas to temperatures of more than 1000 K.

The model procedure is the same as that of Kristensen et al. (2008), except for a change in the geometry. In the original procedure, 1D shock model results were pasted along a parabola to simulate emission from a bow shock observed in H₂ in the Orion Molecular Cloud. Here, in order to simulate C-type shocks along the cavity walls, the 1D shocks are pasted along the ellipsoids from Sect. 2.2.1. The shock velocity is assumed to be constant along the entire length of the wall and the pre-shock density is that of the envelope at the given distance from the protostar. The magnetic field is set to β [2n(H₂)/cm^-3 ] ^1/2, with β fixed at 1 μGauss. Abundances as function of depth into the shock are calculated from a limited chemical network centred on H₂O (Wagner & Graff 1987), which does not include photodissociation or other UV-driven reactions.

The shock model results of Kaufman & Neufeld (1996) are used in combination with the model results of Flower & Pineau des Forêts (2003) and Kristensen et al. (in prep.) to compute the emergent line intensities. The first step is to estimate the extent of the line-emitting region by measuring the FWHM of the cooling-rate profile of the relevant molecule through the shock. The bottom panel of Fig. 3 shows the cooling profiles from Flower and Pineau des Forêts for CO and H₂O, for a shock velocity of 20 km s^-1 and pre-shock densities of 10⁴ and 10^6.5 cm^-3. The corresponding temperature and abundance profiles are plotted in the top panel. CO and H₂O cool the same part of the shock and their cooling lengths are nearly the same. Furthermore, the cooling lengths are nearly inversely proportional to the pre-shock density, i.e., the solid and dotted profiles look similar but are shifted on the logarithmic scale. The cooling lengths provide an estimate of the extent over which the shock is emitting, and as such are an effective measure of the beam-filling factor.

The second step is to divide the cavity wall into discrete segments according to the envelope density. The upper limit of the Kaufman & Neufeld (1996) grid is 10^6.5 cm^-3; the pre-shock density for the denser parts of the envelope is set to that value. Figure 4 shows the adopted segmentation for HH 46; the width of the coloured area is a rough representation of the size of the emitting region. The CO and H₂O line intensities are estimated from the results tabulated by Kaufman & Neufeld. The emission is assumed to fill each of the segments, thus creating an emission map for each line. Just like the raw images produced by LIME (Sect. 2.2.3), the maps are convolved with a Gaussian beam profile of the appropriate diameter. For lines observed with Herschel-PACS, the flux is extracted from the central $\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}\mathrm{\times }\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}$. The fluxes are not corrected for extinction by dust in the envelope.

Our method is specifically designed to simulate optically thin emission from shocks. Low-J CO lines typically observed in and associated with outflows, such as the 2–1 and 3–2, originate in cold (~100 K) swept-up gas present on larger spatial scales than the shocks along the cavity wall (e.g., Arce et al. 2007). These low-J lines are often optically thick and only probe the surface layer of the outflow, where the velocity gradient is largest. Since our focus is on the hot gas observed with Herschel, no attempt is made to quantify the emission from the colder swept-up gas.

An important caveat in the shock models is that they were calculated without a UV field. Several groups are working on irradiated shock models featuring a self-consistent treatment of the UV field in both the pre- and the post-shock gas (e.g., Lesaffre et al. 2011; Kaufman, priv. comm.), but no such models are currently available. Qualitatively, a UV-irradiated shock is expected to be smaller and hotter, and to have a lower H₂O abundance. The latter is a direct effect of the photodissociation of H₂O. The smaller spatial extent and the higher peak temperature are due to the higher degree of ionization in the irradiated pre-shock gas: the ion-neutral coupling length decreases, so the shock front is compressed and the same amount of energy is imparted over a shorter distance. The shock results are discussed in Sect. 4 with these qualitative arguments in mind.

3. Source sample and observations

The three test cases for the model are the low-mass YSOs NGC 1333 IRAS2A, HH 46 IRS and DK Cha (IRAS 12496 − 7650). Each source has been studied extensively in the past, including recent observations with PACS and HIFI on Herschel (van Kempen et al. 2010a,b; Kristensen et al. 2010; Yıldız et al. 2010). The availability of Herschel data and a large body of complementary data is one reason to choose these three sources. Another is that they offer various challenges to the model by spanning a range of evolutionary stages and orientations: IRAS2A is a deeply embedded Stage 0 source seen nearly edge-on (Sandell et al. 1994), HH 46 is a Stage I source seen at a 53° inclination (Nishikawa et al. 2008), and DK Cha is in transition from Stage I to II and is seen nearly pole-on (Knee 1992). The basic observational characteristics and model parameters for each source are summarised in Table 2.

Far-IR spectra of the three sources were obtained with PACS and HIFI on Herschel (Pilbratt et al. 2010; Poglitsch et al. 2010; de Graauw et al. 2010). The HIFI data used in this work were published by Kristensen et al. (2010) and Yıldız et al. (2010). The PACS spectra of HH 46 and DK Cha were published by van Kempen et al. (2010a,b); the PACS spectrum of IRAS2A is previously unpublished. All PACS data were rereduced with the standard pipeline in HIPE v6.1 (Ott 2010). A spectral flatfield was applied to the data to improve the signal-to-noise ratio.

PACS consists of a 5 × 5 integral field unit with $\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}\mathrm{\times }\mathrm{9}\stackrel{\mathrm{″}}{.}\mathrm{4}$ spaxels. For HH 46, the fluxes of the central 9$\stackrel{\mathrm{″}}{.}$4 region are measured directly from the central spaxel and are accurate to 10 − 20%. The observations of DK Cha and IRAS2A were mispointed by a few arcsec. For both sources, the line emission is concentrated near the location of the continuum emission. In order to measure the line fluxes, the spectra in the four spaxels with the brightest continuum emission were first summed. The ratio of the continuum emission in these four spaxels to the total continuum flux in the 5 × 5 array was then used to obtain a wavelength-dependent curve to convert the four-spaxel extraction into a flux in the entire field of view. The fluxes were multiplied by the fraction of the flux that falls in the central spaxel (ranging from 70% below 100 μm to 40% at 200 μm) to estimate the total flux that would have been in the central spaxel of a well-pointed observation. These corrections assume a point source and may therefore overestimate the flux in the central spaxel. Excluding this uncertainty, the fluxes for IRAS2A and DK Cha are accurate to about 30%. Fluxes and upper limits for all three sources are tabulated in Table 3.

The Herschel data are combined with observations of lower-J CO and CO isotopologue lines (up to 7–6) obtained with the James Clerk Maxwell Telescope (JCMT) and the Atacama Pathfinder EXperiment (APEX; Jørgensen et al. 2002; van Kempen et al. 2006, 2009; Yıldız et al. 2010). The different beam sizes for the JCMT, APEX and Herschel are taken into account by convolving the line fluxes from the model with a wavelength-dependent beam size appropriate for the instrument with which each line was observed.

4. Results

4.1. Physical and chemical structure

The best-fit DUSTY envelope model of HH 46 has a radius of 16 100 AU (cut off at 10 K) and a mass of 1.7 M_⊙ (Table 2). The envelope of NGC 1333 IRAS2A is a little smaller (14 100 AU), but because of the shallower density profile its mass comes out 20% higher at 2.0 M_⊙. DK Cha has the smallest and least massive envelope: 5700 AU and 0.02 M_⊙. The density and temperature profiles of each envelope are plotted in Fig. 5.

When the cavities are carved out of the spherical envelopes, the gas and dust are heated up by the stellar UV flux. The best-fit UV luminosities for the three central sources are determined in Sect. 4.2 by fitting a grid of luminosities to the CO observations. The resulting temperature profiles along the cavity walls are shown in Fig. 6 (black and red curves). For comparison, the plot also shows the temperature profiles from DUSTY (blue), which is what the gas and dust would reach in the absence of UV heating. The effect of UV heating is strongest for DK Cha because of the lower envelope densities. The gas temperature peaks at 3800 K at 150 AU from the star. IRAS2A reaches its maximum of 850 K right at the inner edge, followed by a small secondary maximum of 330 K at 320 AU. For HH 46, the density in the inner envelope is high enough that the gas initially remains coupled to the dust and has the same temperature. At larger distances, the decreasing density allows the gas to become hotter than the dust, reaching a maximum of 260 K at 670 AU.

The full 2D gas temperature profiles in Fig. 7 show that the effects of UV heating are not limited to just the cavity wall, but penetrate some depth into the envelope according to the exp( − 0.6A_V) depth dependence from Sect. 2.2.1. The thickness of the UV-heated layer ranges from about 100 AU at r = 1000 AU to about 2000 AU at the outer edge of the envelope. Hence, the UV-heated layer is about a factor of 10 thicker than the shock-heated layer (Fig. 3). Accounting for the predicted effects of the UV field on the shock structure, the thickness ratio between the two layers would be even larger (Sect. 2.3).

In the purely spherical models for IRAS2A and HH 46, the CO abundance follows a drop profile (top panel of Fig. 8 and top left panel of Fig. 9): it is high (~10^-4) in the warm inner parts of the envelope, low (~10^-6–10^-5) at larger radii due to freeze-out onto cold dust, and high again (~10^-4) towards the outer edge because of the low density and corresponding long freeze-out timescale (Schöier et al. 2004; Jørgensen et al. 2005; Yıldız et al. 2010). DK Cha has a somewhat warmer envelope (Fig. 5), preventing CO from freezing out and keeping its abundance at 10^-4 at all radii.

There is no freeze-out either along the cavity wall in any of the 2D models (bottom panel of Fig. 8 and bottom left panel of Fig. 9), because the dust temperature now exceeds the CO evaporation temperature everywhere (Fig. 6). The 10 000 K blackbody adopted for the stellar spectra is energetic enough to dissociate some CO in the cavity wall, reducing its abundance to 10^-6 − 10^-5in IRAS2A and HH 46, and  ~10^-7in DK Cha (Fig. 8, bottom). Moving from the cavity wall towards the midplane, the increasing amount of extinction shields CO from the dissociating photons, and the abundance goes back up to  ~10^-4(Fig. 9; see Fig.  for the other two sources). The dust near the midplanes of HH 46 and IRAS2A remains cold enough for CO to stay frozen, as evidenced by the green region in the bottom left panel of Fig. 9. As a result, horizontal cuts from the cavity wall into the envelope essentially show drop-abundance profiles (Figs. 12 and 13). CO is not expected to freeze out in the less massive and warmer envelope of DK Cha, so the horizontal cuts merely show the effect of the photodissociation rate decreasing with distance to the cavity wall (Figs. 10 and 14).

The protostellar UV field also dissociates H₂O, making for a nearly uniform abundance of 10^-10 along the wall of each of the three sources (Fig. 8). Moving into the envelope, the increasing amount of extinction slows down the photodissociation process and the H₂O abundance increases to a few 10^-8 in the two regions where it does not freeze out: the warm inner envelope and the low-density outer envelope (Fig. 9, bottom right). In between, freeze-out keeps the abundance limited to at most a few 10^-9. As in the case of CO, the original 1D abundance profiles largely survive at the midplane (Figs. 11–14).

With the conditions in the passively heated gas and the UV-heated gas established, the final question is what happens in the small-scale shocks along the cavity walls. As shown in Fig. 3, the peak temperature in a 20 km s^-1 shock ranges from 1300 to 3800 K for pre-shock densities from 10⁴ to 10^6.5 cm^-3. In addition, the peak temperature depends roughly linearly on the shock velocity (e.g., Flower & Pineau Des Forêts 2010). The higher temperatures help overcome the barriers on the O   +   H₂  →   OH   +   H and OH   +   H₂  →   H₂O   +   H reactions (Bergin et al. 1998), so the H₂O abundance in the shocked gas is expected to be orders of magnitude higher than the value of 10^-10 calculated for the UV-heated gas (Fig. 8).

The shock models of Flower & Pineau des Forêts (2003) predict a peak H₂O abundance of  ~10^-4(Fig. 3), but that is in the absence of a UV field (Sect. 2.3). The photodissociation and reformation timescales can be estimated to get an idea of how much the H₂O abundance would change in an irradiated shock. At 1000 AU from the protostar, the UV flux is equivalent to G₀ ≈ 100, which gives a photodissociation timescale of about 1 yr (van Dishoeck et al. 2006). A reformation timescale of about 1 yr occurs for T = 1000 K and n(H₂) = 10⁵ cm^-3 (Bergin et al. 1998). The density at 1000 AU is closer to 10⁶ cm^-3 (Fig. 5), so the reformation timescale is a factor of a few shorter than the photodissociation timescale. However, the two numbers are similar enough that in reality either process could dominate over the other, in particular when considering different points along the cavity wall. A new generation of shock models, incorporating UV radiation, is needed to solve this issue. Until then, the peak H₂O abundance of  ~10^-4in Fig. 3 should be considered an upper limit only.

4.2. UV luminosities and shock velocities

The model has two free parameters – the UV luminosity and the shock velocity – for which the best-fit values are obtained by comparing the computed CO line fluxes to the available Herschel observations in a χ² procedure. The lower-J CO lines observed from the ground are not included in the fit. They contain contributions from the parent molecular cloud and from colder swept-up gas, neither of which is present in our model. The H₂O observations are also not used to constrain the model. The grid of UV luminosities runs from 0.01 to 1 L_⊙. The grid of shock velocities runs from 10 to 35 km s^-1, approximately the critical velocity of C-type shocks for the relevant density regime (Sect. 2.3).

As an example of the dynamic range of the free parameters, Fig. 15 shows the predicted line fluxes for HH 46 for a range of values. The top panel shows how the combined flux from the passively heated envelope and the UV-heated outflow cavity walls changes as function of UV luminosity. The bottom panel does so for the flux from the small-scale shocks as function of shock velocity. These figures are qualitatively the same for the other two sources.

The best-fit UV luminosities and shock velocities are tabulated in Table 2 and shown in Fig. 16. Overplotted in the figure are the 1, 2 and 3σ confidence intervals. Both parameters are well constrained in HH 46. In IRAS2A, the UV luminosity is only constrained from the top. In DK Cha, the shock velocity is unconstrained at the 3σ level and only constrained from below at the 2σ level. Extending the grid of luminosities to lower values does not resolve the issue for IRAS2A, because the gas cannot cool below the dust temperature. The χ² analysis suggests that the CO fluxes observed with Herschel are powered primarily by shocks in IRAS2A, by UV-heating in DK Cha, and by a combination of both in HH 46. This result is discussed in more detail in the next section.

4.3. Line fluxes

4.3.1. Carbon monoxide

With the physical and chemical structure from Sect. 4.1 and the best-fit UV luminosities and shock velocities from Sect. 4.2 (taking L_UV = 0.01 L_⊙ for IRAS2A and ν_s = 35 km s^-1 for DK Cha), the model produces the integrated CO line fluxes shown in Fig. 17. Also plotted are the integrated line fluxes observed from the ground and with Herschel. The only systematic discrepancy between model and observations is an underproduction of the 2–1, 3–2 and 4–3 lines. The same discrepancy was seen by van Kempen et al. (2009) in their purely spherical model of HH 46. The excess flux probably originates in cold gas in the parent molecular cloud or in cold swept-up outflow gas, neither of which is included in our model. The low-J lines of the less abundant C¹⁸O isotopologue are not expected to have any significant cloud or outflow contribution (e.g., Jørgensen et al. 2002; Yıldız et al. 2010), and indeed the model reproduces the available C¹⁸O data to within 50%. The mismatches with some of the individual Herschel data points in Fig. 17 do not appear to follow any systematic trend. They are likely due to uncertainties in the model fluxes and minor factors such as pointing offsets and calibration uncertainties that may be larger than the estimated 10–30%.

The coloured curves in Fig. 17 show the predicted contribution from the passively heated envelope (green), the UV-heated cavity walls (red) and the small-scale shocks (blue) to the total CO emission (black). Each source requires a combination of UV-heated gas and small-scale shocks to reproduce the Herschel data, although there are clear qualitative differences from source to source. IRAS2A is the youngest of the three sources (deeply embedded Stage 0) and its CO ladder appears to be dominated by shock-powered emission. HH 46 is more evolved (Stage I) and shows a 50/50 split between UV-powered and shock-powered emission. DK Cha is still more evolved (in transition from Stage I to II) and its small-scale shocks are responsible for only a small fraction of the overall CO emission. The decreasing contribution of the shock component is consistent with the general trend of the outflow force being weaker in more evolved protostars (Bontemps et al. 1996). The increasing contribution of the UV component is powered by the decreasing densities in the envelope, allowing the gas to become hotter (Figs. 5 and 6; see also Sect. 5.1). Although the sample size is small, the observed data are entirely in keeping with these physically motivated trends. Application of the model to a larger number of sources, such as the full WISH sample, will help to obtain a firmer conclusion. Furthermore, a full parameter study is planned to assess systematically how the CO emission depends on properties such as the envelope mass and the UV luminosity.

The contribution in each curve can be quantified by summing the line fluxes from J_u = 2 to 40 and correcting them for distance to give the overall CO luminosities for the adopted inclinations (Table 4). The CO line fluxes can also be plotted in a rotational diagram to investigate excitation conditions. Table 5 lists the rotational temperatures derived by fitting straight lines to the passive curve from E_u = 10 to 100 K (up to J = 5–4), to the UV curve from 100 to 1000 K (up to J = 18–17), and to the shock curve from 1000 to 4000 K. The error margins reflect the stochastic 1σ uncertainty on each fit, but do not account for any uncertainties in the model results.

The CO ladder for DK Cha in Fig. 17 implies that the six CO lines seen with ISO, from J = 14–13 to 19–18 (Giannini et al. 1999), originate in the UV-heated gas in the outflow cavity walls. Based on a passively heated spherical envelope model, van Kempen et al. (2006) concluded that the ISO lines trace material on larger spatial scales than the material responsible for the 7–6 and lower lines observed with APEX. With the more complex source geometry in our model, such a requirement no longer exists. In fact, the temperature gradient along the cavity wall (Fig. 6), together with the pole-on inclination, suggests the high-J lines seen with ISO originate on smaller spatial scales than the lower-J APEX lines. A more detailed analysis of the spatial extent of the PACS data, using all 5 × 5 spaxels, can help to solve the puzzle, but this is beyond the scope of the current work.

4.3.2. Water

Figure 18 shows the observations and model predictions for H₂O. The complicated rotational structure of H₂O does not allow for plots in the format of the CO ladders in Fig. 17, so the results are shown as rotational diagrams instead. The individual line fluxes can be summed again to obtain the total line luminosities. However, the models only calculate fluxes for a subset of all possible rotational lines. Fluxes for the missing lines in the HIFI and PACS wavelength ranges (up to E_u/k = 2000 K) are calculated by fitting a single rotational temperature to the available lines in each model component (Table 5). The level populations of H₂O are not in LTE, so this method may produce large errors in the flux of a given individual line. However, because those errors tend to cancel when summing over all possible lines, assuming a single temperature is appropriate for calculating the total line luminosities (see also Goldsmith & Langer 1999). The resulting values are tabulated in Table 4.

Unlike CO, the H₂O emission observed with Herschel seems to be dominated (>99%) by shock-heating in all three sources. This is consistent with the broad H₂O line profiles seen with HIFI in IRAS2A and other embedded YSOs (Kristensen et al. 2012). Our shock models actually overproduce many of the detections and upper limits by up to a factor of 10, in particular in HH 46 (Fig. 18). The likely cause is that the H₂O abundances in our shock models were calculated without taking the protostellar UV field into account. Based on the photodissociation and reformation timescales estimated in Sect. 4.1, the H₂O abundance can easily be lowered by a factor of 2 to 10 from the peak value of 10^-4 shown in Fig. 3. At the same time, the abundance of OH would go up, so our shock models are currently expected to underproduce the OH line intensities. However, as argued by van Kempen et al. (2010b) and Wampfler et al. (2010), the OH lines observed with Herschel likely originate in a different physical component than the H₂O lines (a dissociative J-type shock instead of a non-dissociative C-type shock). It is therefore not possible to test whether the OH line intensities predicted for the small-scale shocks are indeed lower than the actual OH emission from that component. A potential problem with lowering the H₂O abundances through photodissociation is that some CO in the shock may be dissociated as well. If that happens, a higher shock velocity is required to fit the CO observations, which in turn would raise the model H₂O fluxes back to their original values of up to an order of magnitude too high. An independent constraint on the CO and H₂O abundances in UV-irradiated shocks is needed to fully resolve this issue.

Another difference between CO and H₂O is seen in the relative contributions of the passive and UV components. For H₂O, both are weak compared to the shock-powered emission, but the passive component is always the brighter of the two (Table 4). The weak emission from the UV-heated gas in the outflow cavity walls is due to the very low H₂O abundance (Figs. 8 and 9), which in turn is due to the wide range of photon energies over which H₂O can be dissociated. CO can only be dissociated by photons of at least 11 eV, so its abundance in the cavity walls remains higher and the UV-heated gas contributes a larger fraction of the overall CO emission.

5. Discussion

5.1. Gas temperature

The largest unknown in our model is the gas temperature due to UV heating of the outflow cavity walls. Calculating the gas and dust temperatures for a given density and UV flux is a standard problem in models of photon-dominated regions (PDRs; e.g., Tielens & Hollenbach 1985; Hollenbach & Tielens 1997). However, the PDR community has not yet converged on a unique solution: in a benchmark test of ten different PDR codes, the temperature solutions varied by up to an order of magnitude (Röllig et al. 2007). As shown in Sect. 4.3, the temperature grid from Kaufman et al. (1999, hereafter K99), along with the exp( − 0.6A_V) depth dependence derived from the benchmark test, works well to reproduce the CO observations. However, given the complicated nature of our model, this does not necessarily mean that the K99 grid gives the correct temperatures. For example, if the grid systematically overpredicts the temperature, it may still give a good match to the data by decreasing the UV luminosity. Other temperature grids may also be able to match the data. The goal of this section is two-fold: (1) to show that, at least on a qualitative level, the K99 grid follows the density and UV flux dependence expected from basic heating and cooling physics for the relevant range of conditions; and (2) to see how robust our results are to changes in the adopted temperature grid.

The K99 grid is reproduced in Fig. 19. In its original form, it only covers H₂ densities up to 10⁷ cm^-3. The densities in our model go up to 10⁹ cm^-3, so the grid is extrapolated as T_gas ∝ 1/n(H₂). The n^-1 dependence comes from the approximate functional forms of the heating and cooling mechanisms expected to dominate at high densities: photoelectric heating ( ∝ n; Bakes & Tielens 1994) and gas-grain collisional cooling ( ∝ n²; Hollenbach & McKee 1989). The density dependence of the temperature in the original K99 grid is almost exactly n^-1 between 10⁶ and 10⁷ cm^-3.

Overplotted as black lines in Fig. 19 are the gas density and incident UV flux along the cavity walls of the three sources. The density and UV flux both fall off approximately as r^-2, so they scale linearly with each other. At densities above  ~10⁶ cm^-3, the gas temperature is expected to scale with the inverse of the density (see above) and linearly with the UV flux (Bakes & Tielens 1994), and should therefore remain constant along the inner part of the wall. This is indeed seen in Fig. 6: T_gas does not change much over the range of radii corresponding to densities of more than 10⁶ cm^-3. Photoelectric heating likely remains dominant at lower densities, but atomic fine-structure lines take over from gas-grain collisions as the dominant cooling mechanism. Their cooling rate depends linearly on the density as long as the density exceeds the critical density of 3 × 10³ cm^-3 for the [C ii] 158 μm line, so the temperature is expected to remain roughly constant with density (K99; Meijerink et al. 2007). Because the G₀ dependence remains, the temperature along the cavity wall should decrease as r^-1 for n(H₂) < 10⁶ cm^-3. Again, this is the approximate behaviour seen in Fig. 6. The dependencies change again for densities below 3 × 10³ cm^-3, but such low densities do not appear in our model.

The location of the black lines in Fig. 19 shows why the UV-powered CO emission looks qualitatively different in DK Cha than in IRAS2A and HH 46. The lower envelope densities in the former turn the conditions along the cavity wall into a classic PDR, with densities of 10⁶–10⁷ cm^-3 and UV fluxes of 10³ − 10⁴ times the mean interstellar flux. The gas therefore gets much hotter and CO is excited to higher rotational levels. The steep temperature gradient in this part of Fig. 19 also means that the CO emission is rather sensitive to details in the SED fitting with DUSTY. A factor of two difference in density leads to a similar change in the gas temperature, which can affect the highest-J CO lines by a few orders of magnitude. Figure 17 tentatively shows an evolutionary trend of the CO emission in more evolved sources being more strongly UV-powered. While the application of a physical model is necessary to quantify exactly how much of the total CO ladder in a given source is due to UV heating, the observations are consistent with the main point that with decreasing density (assuming sufficient UV), sources becomes more PDR-like. This is what one would naturally expect in the transition from Stage 0 to Stages I and II objects as presented in this paper.

Of equal importance to the CO emission is the choice of the PDR temperature grid. An easy way to get a rough idea of how changes in the adopted grid would affect the CO line fluxes is to parameterise the gas temperature based on the dominant heating and cooling mechanisms. Two equations suffice to do so. The first one approximates the gas-grain collisional cooling rate (Hollenbach & McKee 1989): $\begin{array}{ccc}{\mathrm{\Lambda }}_{\mathrm{g}\mathrm{-}\mathrm{g}}& \mathrm{=}& \mathrm{(}\mathrm{1}\mathrm{\times }{\mathrm{10}}^{-31}\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathrm{)}{\mathit{n}}_{\mathrm{H}}^{\mathrm{2}}\sqrt{\frac{{\mathit{T}}_{\mathrm{gas}}}{\mathrm{1000}\hspace{0.17em} \mathrm{K}}}\sqrt{\frac{\mathrm{100}\hspace{0.17em}\AA }{{\mathit{a}}_{\min }}}\\ & & \mathrm{\times }\left[\mathrm{1}\mathrm{-}\mathrm{0.8}\exp \left(\mathrm{-}\frac{\mathrm{75}\hspace{0.17em} \mathrm{K}}{{\mathit{T}}_{\mathrm{gas}}}\right)\right]\mathrm{(}{\mathit{T}}_{\mathrm{gas}}\mathrm{-}{\mathit{T}}_{\mathrm{dust}}\mathrm{)}\mathit{,}\end{array}$(6)where n_H is the total hydrogen nucleus density (equal to 2n(H₂) in our model for all practical purposes) and a_min is the minimum grain size (set to 10 Å). The second equation is for the photoelectric heating rate (Bakes & Tielens 1994): ${\mathrm{\Gamma }}_{\mathrm{pe}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{\times }{\mathrm{10}}^{-24}\hspace{0.17em}\mathrm{erg}\hspace{0.17em}{\mathrm{cm}}^{-3}\hspace{0.17em}{\mathrm{s}}^{-1}\mathrm{)}\mathit{ϵ}{\mathit{G}}_{\mathrm{0}}{\mathit{n}}_{\mathrm{H}}\mathit{.}$(7)The photoelectric heating efficiency, ϵ, is treated as a free parameter for the purpose of this simple approximation. As such, it hides any uncertainties in the grain charge and in the numerical prefactors in the two equations (see, e.g., Young et al. 2004). Given a dust temperature from RADMC (or, alternatively, from the analytical expression of Hollenbach et al. 1991), the gas temperature at n(H₂) > 10⁶ cm^-3 is obtained by solving for Λ_g − g = Γ_pe. At lower densities, the dependence of T_gas on n(H₂) disappears to zeroth order, and we set T_gas [G₀,n(H₂) < 10⁶  cm³ ]  = T_gas [G₀,10⁶  cm³ ] .

Figure 19 shows the 10, 30, 100 and 1000 K contours for ϵ = 0.2 as dashed white lines. The approximate grid reproduces the K99 grid to within a factor of 2 for the full range of densities and UV fluxes encountered along the three cavity walls. The approximation breaks down for G₀/n(H₂) < 10^-4 (high density, low UV) and G₀/n(H₂) > 10^-3 (low density, high UV), and cannot be used in either of those regimes. The value of 0.2 required for a good match with K99 exceeds the maximum value of 0.05 allowed by Bakes & Tielens (1994). However, the approximate temperature grid is purely meant as an easy tool to study the effect of different temperature grids, so this factor of 4 should not be considered physically important. The key point is that the K99 grid and the approximate grid show the same qualitative trends for the relevant range of conditions, which means the K99 grid follows the behaviour expected from the dominant heating and cooling mechanisms.

With ϵ = 0.2, the approximate grid does a good job of reproducing the CO line emission for the UV component. Shown in Fig. 20 are the integrated line fluxes from the UV-heated gas in HH 46, using either the K99 grid (dashed black curve, identical to the red curve from Fig. 17) or the approximate grid with three different values of ϵ (dashed coloured curves). The ϵ = 0.2 curve lies within 30% of the K99 curve up to J_u = 20, at which point shocks take over as the dominant producer of CO emission. When adding the passive and shock components to the ϵ = 0.2 curve (producing the solid red curve), the overall agreement with the data is as good as in Fig. 17. Taking a lower or higher photoelectric heating efficiency would under- or overproduce the J = 14–13, 16–15 and 18–17 lines.

The photoelectric heating rate scales linearly with ϵ and G₀ (Eq. (7)), so the three values of ϵ in Fig. 20 can be brought to give the same CO fluxes by changing the UV luminosity. Herein lies the degeneracy alluded to at the beginning of this section and in Sect. 4.2: as long as two temperature grids show the same qualitative trends, both can reproduce the CO observations, but they have different best-fit UV luminosities. This degeneracy can only be lifted by calibrating the temperature grids against a source whose incident UV flux can be constrained independently in the same part of n(H₂)–G₀ parameter space.

5.2. Caveats in the shock models

The shock modelling relies on three assumptions: (1) the shocks move along the cavity walls; (2) the effective beam-filling factor can be obtained from accurately estimating the shock width; and (3) the protostellar UV field does not affect the H₂O abundance. In order to test where the emission is coming from in this scenario – extended lower-density gas with a large beam-filling factor or higher-density gas with a small beam-filling factor – fluxes from the standard shock model are compared to fluxes from models without any contribution from gas with densities of 10⁴ cm^-3 and lower, 10⁵ cm^-3 and lower, or 10⁶ cm^-3 and lower. The shock velocity is kept at 20 km s^-1 for this experiment. For the case of HH 46, gas at densities below 10⁴ cm^-3 is found to contribute less than 10^-4%, and gas at densities between 10⁴ and 10⁵ cm^-3 only 0.3%. About 75% of the overall flux comes from the segment with a pre-shock density of 10^6.5 cm^-3. The other two sources show almost the same numbers, indicating that the shocked gas observed by PACS is dominated by emission from the densest parts of the YSO.

If the density is assumed to be fixed at more than 10⁶ cm^-3, the shock velocity alone determines the absolute line fluxes. This was already seen in the bottom panel of Fig. 15, where the basic shape of the CO ladder is nearly independent of velocity. Hence, the velocity determines the intensity of the CO emission, and the pre-shock density determines the shape of the ladder. Comparing the fluxes amongst the three sources, the differences are dominated by their respective distances. This confirms that the actual geometry only plays a small role with respect to the predicted CO emission, and that the bulk of emission is created in very dense gas contained within the central PACS spaxel.

If the effective beam-filling factor is treated as a free parameter, it becomes degenerate with the shock velocity. A smaller beam-filling factor would require a higher shock velocity to produce the same line fluxes, and vice versa. However, for a given source geometry, such as the ellipsoid outflow cavity in our model, the beam-filling factor follows from the cooling lengths and is not truly a free parameter. The only way to constrain either the geometry or the beam-filling factor would be to have an independent handle on the pre-shock density and the shock velocity. Our method therefore gives a unique solution for the given shock geometry, but other geometries cannot be excluded.

6. Conclusions

The Herschel Space Observatory offers a new probe into the hot gas present in embedded low-mass young stellar objects (YSOs) through emission lines of highly rotationally excited CO (up to 4000 K above the ground state) and H₂O (up to 600 K). This paper presents the first model to reproduce these observations quantitatively. The model starts with a passively heated spherical envelope, which fits the lowest-J CO isotopologue lines observed from the ground, but falls short of the ${{12}}^{}$CO Herschel data by several orders of magnitude. Two additional model components are therefore invoked: (1) gas in the walls of a bipolar outflow cavity, exposed to the protostellar UV field; and (2) small-scale C-type shocks along the cavity walls.

The main conclusion is that, within the context of the model, the far-infrared CO rotational line emission in the low-mass YSOs NGC 1333 IRAS2A, HH 46 and DK Cha can only be reproduced with a combination of UV-heated gas and C-type shocks. The three sources show a tentative evolutionary trend: the CO emission appears to be dominated by shock-heating in the youngest source (IRAS2A), by UV-heating in the oldest source (DK Cha), and by a 50/50 combination in the source of intermediate age (HH 46). This trend is powered by the envelope density decreasing as a protostar evolves from Stage 0 to Stages I and II.

The H₂O lines are dominated by shocks in all three sources, with less than 1% of the total H₂O line luminosity originating in UV- and passively heated gas. This is a direct result of the different abundances in the three components. Freeze-out keeps the H₂O abundance below 10^-8 in the cold outer envelope. The dust grains in the UV-heated layer are warm enough for H₂O to evaporate, but it is dissociated by the UV field once it does. Efficient reformation at high temperatures boosts the abundance to 10^-5−10^-4in the C-type shocks, so they are responsible for most of the H₂O emission. CO has a lower evaporation temperature (20 versus 100 K), so its abundance is orders of magnitude higher than that of H₂O in the passively heated and the UV-heated gas, allowing these two components to contribute more strongly to the overall CO emission.

The main uncertainty in the model lies in the calculation of the gas temperature as function of the incident UV flux. This is a well-known problem in photon-dominated regions (PDRs), and many solutions are available in the literature. However, they disagree by up to an order of magnitude and some do not even agree qualitatively on how the gas temperature depends on the density or the UV flux. The sensitivity of our model results to the treatment of the gas temperature is tested with an approximate analytical description of the expected dominant heating and cooling mechanisms. This exercise shows that the quantitative contributions of the UV-heated gas and the C-type shocks to the hot gas emission observed by Herschel are uncertain by at least a factor of 2, but the qualitative conclusion that both components are needed to explain the observations is robust.

Online material

---

Acknowledgments

The authors are grateful to the WISH and DIGIT teams – in particular Neal Evans, Jes Jørgensen and Ted Bergin – for stimulating discussions and scientific support. Constructive comments from the anonymous referee are kindly acknowledged. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by a Spinoza grant and grant 614.001.008 from the Netherlands Organisation for Scientific Research (NWO), and by the European Community’s Seventh Framework Programme FP7/2007–2013 under grant agreement 238258 (LASSIE). Partial support for this work was provided by NASA through RSA awards No. 1371476 and 1416374, by the NSF through grant AST-1008800, and by JPL through Contract 1358118. HIFI has been designed and built by a consortium of institutes and university departments from across Europe, Canada and the United States under the leadership of SRON Netherlands Institute for Space Research, Groningen, the Netherlands, and with major contributions from Germany, France and the USA. Consortium members are: Canada: CSA, U. Waterloo; France: CESR, LAB, LERMA, IRAM; Germany: KOSMA, MPIfR, MPS; Ireland: NUI Maynooth; Italy: ASI, IFSI-INAF, Osservatorio Astrofisico di Arcetri-INAF; Netherlands: SRON, TUD; Poland: CAMK, CBK; Spain: Observatorio Astronómico Nacional (IGN), Centro de Astrobiología (CSIC-INTA); Sweden: Chalmers University of Technology (MC2, RSS & GARD), Onsala Space Observatory, Swedish National Space Board, Stockholm University (Stockholm Observatory); Switzerland: ETH Zurich, FHNW; USA: Caltech, JPL, NHSC. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain).

References

All Tables

All Figures