© ESO, 2013

1. Introduction

Because young high-mass stars (≳8 M_⊙) reach the hydrogen-burning phase before the accretion phase ends (Palla & Stahler 1993), their formation process differs from less massive protostellar objects. Several theories exist for the formation of massive stars (see review by Zinnecker & Yorke 2007), including stellar mergers and so-called “competitive accretion” (Bonnell et al. 2001; Bonnell & Bate 2005). However, the accumulating observational evidence for collimated outflows and disk structures in massive star-forming regions (Beuther et al. 2002; Cesaroni 2005; Sandell & Wright 2010; Kraus et al. 2010; Wang et al. 2012) supports a picture of anisotropic accretion through a disk structure, analogous to the low-mass star formation process (e.g., Shu et al. 1987; Larson 2003), at least for luminosities up to 10⁵ L_⊙. Even with the latter view of high-mass star formation, striking differences arise compared to low-mass star formation theory: shorter timescales by a factor ~10 and higher internal luminosities by several orders of magnitude. Moreover, the bulk of the energy produced by massive stars is emitted at UV wavelengths, by which the surrounding gaseous material is heated, partly ionized, and otherwise chemically influenced. The effect of UV radiation on circumstellar gas has been recognized for low-mass protostars (e.g., Spaans et al. 1995) and is even stronger in the high-mass case (Stäuber et al. 2004). Since a collapsing protostellar envelope depends on gas-phase molecules and atoms to dissipate gravitational energy, knowledge of the chemical composition and excitation conditions of atomic and molecular gas is key in the understanding of the high-mass star formation process.

AFGL 2591 is a high-mass star-forming region with a centrally heated molecular envelope. It is located at galactic coordinates in the Cygnus-X region, which has an average distance of 1.7 kpc (Motte et al. 2007). While the distance to AFGL 2591 itself has been uncertain in the past (Van der Tak et al. 1999; Schneider et al. 2006), in this paper we adopt a distance of 3.3 kpc as derived from trigonometric parallax measurements of H₂O masers by Rygl et al. (2012). With respect to the distance of 1 kpc often used in previous work (e.g., Hasegawa & Mitchell 1995; Van der Tak et al. 2000; Van der Tak & Menten 2005; Van der Wiel et al. 2011), the new distance changes the estimates for the mass of the dominant protostar, VLA 3, to ~40 M_⊙ and the bolometric luminosity to 2 × 10⁵ L_⊙ (Sanna et al. 2012; Jiménez-Serra et al. 2012). All linear scales mentioned in the remainder of this paper are based on a distance of 3.3 kpc to AFGL 2591, as is the numerical model description we use in Sect. 4. AFGL 2591 exhibits a powerful bipolar molecular outflow extending to > 1′ from the central source (Lada et al. 1984; Hasegawa & Mitchell 1995). At a 3.3 kpc distance, this corresponds to a linear scale ~1 pc, which is not untypical for massive protostellar outflows (Beuther et al. 2002). A massive disk has been proposed to exist around source VLA 3, but with an angular extent of < 1′′ (Van der Tak et al. 2006; Wang et al. 2012) the disk will not be resolved in the current study, which traces ≳ 10′′ scales. Although interferometric observations have revealed multiple sources inside the envelope at spacings of ≲ 10 000 AU (Van der Tak et al. 1999; Trinidad et al. 2003; Benz et al. 2007; De Wit et al. 2009; Sanna et al. 2012), AFGL 2591 is isolated at scales of a few parsecs. Such isolation is exceptional for massive star-forming envelopes, which are usually found in clustered environments. The absence of complex dynamical and radiative interaction with neighboring objects has made AFGL 2591 a popular high-mass protostellar object to study.

Most previous modeling efforts for the large scale envelope of AFGL 2591 have assumed a spherical morphology with a power law density profile (Van der Tak et al. 1999, 2000; Doty et al. 2002; Stäuber et al. 2005; Benz et al. 2007; De Wit et al. 2009). In these studies, gas temperatures range from ~20 K in the outskirts of the envelope (~10⁵ AU from the center) to several hundred K at distances as close as a few hundred AU from the protostar. Likewise, Van der Wiel et al. (2011) employ spherically symmetric model geometries, but they also turn to flattened axially symmetric geometries to explain molecular emission maps at scales of ~10⁴–10⁵ AU. Moreover, after suggestions for the presence of cavities in the envelope carved out by outflow motions (Van der Tak et al. 1999), and observational evidence in the near-infrared by Preibisch et al. (2003), detailed radiative transfer models including outflow cavities are presented by Bruderer et al. (2009a,b, 2010a). They include a chemical balance to allow for photodissociation of molecules along the interface between the envelope and the cavity, and show that a considerable fraction of the gas is at temperatures exceeding 1000 K. This adjusted temperature balance and the possible presence of shocked regions obviously affect the excitation conditions of the molecular gas.

In the envelopes of massive star-forming regions, much of the gas is well below 50 K, but temperatures also exceed 100 K in a significant portion of the gas. In this work, we study this warm gas in the large scale envelope of AFGL 2591. To study these high gas temperatures with ground-based observations, one relies on low abundance species with relatively complex spectra such as CH₃OH, SO₂ and H₂CO. Because H₂O is only abundant at temperatures above ~100 K (Ceccarelli et al. 1996; Doty & Neufeld 1997), the most important probes of star-forming gas are linear rotor molecules. These species, however, are only observable from the ground in transitions that trace either very hot (~10³ K for near-infrared H₂ lines) or cold gas (≲50 K for low J transitions of asymmetric linear rotors in the (sub)millimeter). Space-based observations are therefore necessary to obtain a comprehensive picture of the warm gas at a few hundred K. We exploit the far-infrared frequency range accessible to the Herschel Space Observatory (Pilbratt et al. 2010) and its Heterodyne Instrument for the Far-Infrared (HIFI, De Graauw et al. 2010) to observe CO, HCO⁺, CS, HCN and HNC in high rotational transitions with upper level energies ranging from 80 to ~700 K.

This paper is structured as follows. Section 2 describes the observations and processing of spectra obtained with Herschel/HIFI from space and of complementary ground-based millimeter-wave observations. In Sect. 3 we present line detections of CO and isotopologues, HCO⁺, CS, HCN and HNC, and we discuss trends in line position, shape and intensity. The observed line strengths of the envelope component are compared with radiative transfer models in Sect. 4. Section 5 presents a derivation of physical conditions in kinematically distinct components. Finally, we discuss the results in Sect. 6 and summarize our conclusions in Sect. 7.

2. Observations

2.1. 490–1850 GHz Herschel/HIFI spectroscopy

The Herschel Key Program “Chemical HErschel Surveys of Star-forming regions” (CHESS, Ceccarelli et al. 2010) uses the HIFI instrument to perform spectral surveys of various Galactic star-forming regions. One of the targets is AFGL2591, for which an unbiased spectral survey was conducted in the 490–1240 GHz range (bands 1a−5a, totaling 18.4 h of observing time), complemented by selected frequency settings between 1260 and 1850 GHz (bands 5b–7b). Of the total of eight targeted frequency settings, this work only uses the five that cover ¹²CO transitions (2.1 h of observing time).

The observed frequency bands are listed in Table 1, along with observing dates, typical values for the noise level, system temperature, and half-power beam width (HPBW). Observations were carried out using the dual beam switch mode with off positions 3′ east and west of the target position: RA = 20^h29^m249, δ = 40°11′21′′ (J2000). The spectral scans (bands 1a−5a) were conducted with the wideband spectrometer (WBS) using a redundancy factor of 4. Per band, per polarization (H or V), per sideband (lower or upper), the scans contain 90–151 datasets, each composed of four 1 GHz “subband” spectra. The pointed frequency settings (bands 5b–7b) were also executed in the dual beam switch mode, in this case with the fast chop and stability optimization options selected. The spectral resolution of the WBS is 1.1 MHz, corresponding to 0.66 km s^-1 at 500 GHz and 0.18 km s^-1 at 1850 GHz.

The raw data of the spectral scans were processed with the standard pipeline included in HIPE (Ott 2010) version 8.1. Further processing of the individual datasets in the spectral scans was done in HIPE, using scripts written by the CHESS team (Kama et al. 2013): (1) data inspection and flagging of spectral regions with spurious features; (2) subtraction of sinusoid standing waves; (3) subtraction of a polynomial baseline, generally of order ≤ 3, fitted to line-free regions. Strong spectral lines are known to create minor ghost lines in the sideband deconvolution process (see below). Therefore, the above steps were repeated by flagging not just the spurious features, but also the strong spectral lines, resulting in a cleaner secondary spectrum with all strong lines masked. The threshold for marking lines as strong varied from T_A^∗ = 1 K in band 1 to 8 K in band 5a, depending on the amount of lines and the noise level in a particular band.

The spectral scans of the CHESS program were executed such that each sky frequency was covered by several settings of the local oscillator, which allows for the mathematical disentanglement of the upper and lower sideband signals that result from double sideband heterodyne observations (Comito & Schilke 2002). This “sideband deconvolution” was performed with sideband gain ratios fixed to unity. Spectra from the HIFI spectral scans (490–1240 GHz) presented in this paper are all single sideband, “unfolded” spectra.

The correspondence of deconvolved products for the two polarization backends, H and V, is checked and found to agree to within 3% in peak intensity for various lines, and line peak positions and line shapes generally match to within the noise level (cf. Table 1). The signal from the two polarizations is averaged by feeding the datasets from the two backends to the sideband deconvolution algorithm simultaneously. One exception is band 4b, where we used only the H polarization, since V is contaminated by numerous spurious features and the V-only deconvolved spectrum exhibits noise levels of up to four times higher than in the H polarization.

The WBS spectra from the pointed frequency settings of CO transitions above 1240 GHz (bands 5b, 6a, 6b, 7a and 7b) are processed with HIPE version 7. Subsequent removal of polynomial baselines of order ≤ 4 and sine-shaped standing waves is applied in the CLASS package¹. With a single frequency setting, sideband deconvolution is impossible, so the high-frequency CO lines are double sideband spectra. Like in most spectral scans described above, the average of the two polarizations is used, with the exception of CO 11–10 in band 5b. For this line, the spectral line shape shows discrepancies between H and V. Since the baseline rms noise level is ~20% higher in V, we only use the H polarization in this case.

We use the main beam efficiency η_mb, as determined from a beam characterization experiment based on Mars observations (Roelfsema et al. 2012), to convert the T_A′ intensity scale to main beam intensity: T_mb = (η_fwd/η_mb)T_A^∗ with the forward efficiency η_fwd = 0.96 for HIFI. Although η_mb is in principle frequency-dependent, it varies by < 1% across an individual band. A more significant variation in the value of η_mb is the 14% drop between band 4b and 5a (Table 1).

The scaling to main beam intensity is strictly only applicable to emission regions that are not larger than the main beam of the telescope. For the lines studied here, this criterion is largely obeyed, so we apply the η_mb factor to facilitate a direct comparison to modeled intensities.

Conservative estimates put the formal uncertainty for the intensity calibration of HIFI at 9% in bands 1 and 2, at 10% in bands 3 and 4, and ~14% in bands 5–7 (Roelfsema et al. 2012). However, the uncertainty budget is often dominated by standing waves and the uncertainty in the sideband ratio, both of which are mitigated considerably by our systematic cleaning (all bands) and sideband deconvolution of the spectral scans (bands 1a−5a). Based on consistency checks of HIFI line measurements over Herschel’s operational lifetime, the observed depth of saturated absorption lines, comparison to ground-based observations and cross-comparison with Herschel/SPIRE, realistic values for HIFI’s absolute calibration uncertainty should be between 5 and 10%. We adopt a uniform uncertainty of 10% for the line intensities measured from our HIFI spectra.

2.2. Ground-based (sub)-millimeter spectroscopy

Apart from the Herschel spectra, this paper uses ground-based observations with the IRAM 30 m Telescope² and the James Clerk Maxwell Telescope³ (JCMT).

IRAM 30 m observations were performed on July 21, 2010. The observations were conducted with the EMIR0 and EMIR2 receivers in parallel (two linear polarizations each) with the VESPA autocorrelator backend in the 2.7 mm (~110 GHz) and 1.3 mm (~230 GHz) windows⁴. The C¹⁸O and ¹³CO J = 1−0 lines were observed simultaneously by splitting the backend in two subwindows centered on the lines with a frequency sampling of 78 kHz (~0.2 km s^-1). The C¹⁸O 2–1 line was observed with the same frequency sampling (~0.1 km s^-1). We used frequency switching mode and an integration time of 1 min on the sky (2 min after averaging the two polarizations). The weather was good with precipitable water vapor of 3 mm, and system temperature was 100 K at 2.7 mm and 200 K at 1.3 mm, in T^∗_A scale. A fourth order polynomial baseline has been subtracted. The measured noise in the spectra, in T_mb units, is 46 mK for the 2.7 mm spectra and 116 mK for the 1.3 mm spectrum after averaging and folding. We use HPBW of 22′′ and 11′′, and main beam efficiencies (η_mb) of 78% and 61% for 2.7 mm and 1.3 mm, respectively. For the IRAM/EMIR spectra, we take an absolute intensity calibration uncertainty of 10%.

Spectra of rotational transitions of CO species, HCN, HNC, HCO⁺, and CS that fall between 330 and 373 GHz are extracted from JCMT Spectral Legacy Survey data (Plume et al. 2007; Van der Wiel et al. 2011). The beam size of the JCMT is 14–15′′ at these frequencies, and the spectral resolution used in the survey is 1 MHz, corresponding to ~0.8 km s^-1. The absolute intensity calibration for the JCMT measurements is uncertain to 15%. A detailed description of the data acquisition and processing is given by Van der Wiel et al. (2011).

3. Observational results

3.1. Detected emission lines and line shapes

The spectral scans and targeted frequency settings described in Sect. 2 cover high-J transitions of the linear molecules HCO⁺, CS, HCN, HNC, CO, ¹³CO, C¹⁸O, and C¹⁷O, as listed in Table 2. Other detected spectral lines will be presented in a companion paper by Kaźmierczak et al. (in prep., hereafter “Paper II”).

For all but one of the species mentioned above, frequencies of rotational transitions are extracted from the Cologne Database for Molecular Spectroscopy (CDMS, Müller et al. 2005). For C¹⁷O we use the JPL frequencies (Pickett et al. 1998)⁵. We search the processed HIFI spectra for line signal near these frequencies. For weak lines (peak intensity <5σ_rms) the spectrum is smoothed in steps of a factor 2, until channels become so wide that they become comparable to the line width. If the S/N is still below 5 at this point, a conservative upper limit to the integrated line flux is then set at 5σ_rms × (3 km s^-1), with the last factor motivated by the typical line width expected from firmly detected lines at lower J.

Spectral line profiles of transitions of CO, ¹³CO, C¹⁸O, and C¹⁷O are presented in Fig. 1; those of HCO⁺, CS, HCN, and HNC are shown in Fig. 2 (see also Table 2). If a line needs smoothing of the spectrum before its peak intensity is above 5σ_rms, Fig. 2 gives the smoothed spectrum.

All lines with sufficient signal-to-noise show an asymmetric profile with line wings that are stronger on the blueshifted than on the redshifted side. The blueshifted line wing (V_lsr ≲ −10 km s^-1) is visible in emission in the lowest HIFI transitions of CO, and to a lesser degree of HCO⁺. We attribute the asymmetry on the blue side of the line profile to the approaching outflow (see Sect. 3.2). In both CO and HCO⁺, it becomes less pronounced at higher J-levels, but it is nevertheless present up to the highest CO transition observed (J = 16−15). A separation between envelope and outflow emission is also seen in CO 2−1 spectra by Mitchell et al. (1992), similar to our CO 3–2 profile. In addition, this line wing is known from infrared absorption studies of ¹²CO and ¹³CO. In ¹²CO, the outflow is seen out to 200 km s^-1 blueshifted (Mitchell et al. 1989; Van der Tak et al. 1999).

Fig. 1 Spectral line profiles of CO, 13CO, C18O, and C17O. The 1–0 and 2–1 spectra are from the IRAM 30 m telescope, the 3–2 lines are from JCMT/HARP-B. All higher-J lines are from Herschel/HIFI. For clarity, consecutive spectra, each labeled by their upper and lower J-level, are offset in the vertical direction. Corresponding rest frequencies and energy levels are listed in Table 2. The systemic velocity ( − 5.8 km s-1) of AFGL 2591 is indicated by the dashed vertical line.

Fig. 2 As Fig. 1, but for HCO+, CS, HCN, and HNC. All lines are measured with Herschel/HIFI, except HCO+, HCN, and HNC 4–3, as well as CS 7–6, which are measured with JCMT.

Contrary to most bright lines, which peak at the systemic velocity and have a blueshifted line wing, the peaks of all ¹²C¹⁶O lines (Fig. 1) are displaced toward redshifted velocities. This is attributed to the optical depth and the possible presence of absorption components affecting the line shape. Both these effects scale with the overall abundance, which is far higher for ¹²C¹⁶O than for any other species studied here. In addition, the 12–11 and 14–13 lines of CS (Fig. 2) appear asymmetric with excess emission on the redshifted side at the ~2-σ level.

The 11–10 line of HCO⁺ appears to be split in two separate components. The brightest one follows the trend of increasingly redder peak velocity with higher energies, while a secondary, narrow peak emerges at −8 km s^-1 (5σ detection). The same feature, although fainter than for J = 11–10, is visible in the 10–9 and possibly also the 9 − 8 lines of HCO⁺ when a single component profile based on lower lines is subtracted. In the other species studied here, this additional narrow component is found only in C¹⁷O 8–7 at the 2.5-σ level, while it might be hidden in the line wing of C¹⁷O 7–6. Based on the firm detection in HCO⁺, we conclude that the feature must originate from dense and very warm gas, given that it is brighter at J = 11–10 than at lower transitions.

3.2. Observed spectral line energy distributions

Integrated line intensities, , are determined by fitting Gaussian profiles to the observed spectral lines (Table 2), using the Levenberg-Marquardt fitter in the line analysis module of cassis⁶. First, a single Gaussian profile is fitted to each observed line⁷. In cases where this yields a residual spectrum with an rms noise comparable to or below the local baseline noise level, a single Gaussian profile is adopted. If subtraction of the single Gaussian fit leaves a significant residual signal, a fit of two blended Gaussian profiles is performed. If this gives a residual spectrum with a significantly lower rms than the one resulting from the single Gaussian fit, we adopt the double Gaussian fit. Examples of cases where a single or a double Gaussian fit is appropriate are shown in Fig. 3. Although the CO lines up to J = 5–4 may accommodate more than two Gaussian components, the additional components are so faint that they have negligible effect on the two primary components.

Fig. 3 Examples of Gaussian component fits to HCN 6–5 (top, single component) and HCO+ 6–5 (bottom, two components). Observed line profiles are shown as histograms, fit components are shown as smooth solid curves and the sum of the two components in the bottom panel is represented by the dashed curve.

Lines that are well represented by a double Gaussian fit are CO, ¹³CO, C¹⁸O and HCO⁺ up to a certain J level, as well as all ground-based line profiles, while a single Gaussian suffices for the HIFI lines of C¹⁷O, HCN, HNC, and CS. In general, a lower signal-to-noise of a spectral line reduces the chance of finding a meaningful secondary Gaussian component. For the lines with a double Gaussian profile, we attribute the component with the narrowest line shape and reddest V_centr to the quiescent envelope, and the broader blueshifted component to outflow material. For lines with a single Gaussian profile, all emission is attributed to the envelope. For CO and HCO⁺, the component separation allows for an analysis of variations in physical conditions between the envelope and the outflow. This is addressed in Sect. 5.2.

Figure 4 shows the of the envelope components versus the energy of the upper level (E_up). For the lines where only a single Gaussian fit is used, it is still likely that up to ~50% of this emission stems from the outflow component. We use a gray errorbar to indicate the intensity level that remains after 0–50% of the single Gaussian intensity is removed from the envelope component intensity. The general trend in these spectral line energy distributions (SLEDs) is that the line intensities decrease with increasing upper level energy. A detailed comparison of the observed envelope SLEDs to models is discussed in Sect. 4.

We emphasize that, although we present various trends as a function of E_up, it is important to recognize the additional systematic effects of decreasing beam size, and decreasing line optical depth, simultaneously with increasing E_up for consecutive transitions.

3.3. Observed centroid velocities and line widths

From the fitted Gaussian profiles, we also obtain measures for centroid velocity and line width (Table 2). Figure 5 shows the trend of fitted values of V_centr and FWHM for the envelope component, as a function of the upper level energy of the line transition. We only plot those species for which at least two lines are fitted with a double Gaussian profile. The main form of CO is excluded, since its line widths are heavily influenced by optical depth effects. Results from the HIFI spectral scans of the low-mass protostar IRAS 16293, also from the CHESS program, show a similar trend as is found in our data of AFGL 2591, but with V_centr decreasing with increasing E_up (Emmanuel Caux, priv. comm.). We therefore conclude that the shifting centroid velocity over the HIFI bands is not an instrumental effect.

In the diagram for V_centr (Fig. 5, top), the two CO isotopologues appear to have a minimum V_centr (bluest peak) at E_up ~ 80 K (J = 5–4), with a monotonously increasing trend toward higher E_up. The four HCO⁺ lines also exhibit an overall trend of centroid velocities becoming redder with increasing E_up. To explain this systematic shift in centroid velocity, we note that increasingly higher J levels do not only probe higher temperatures, parameterized by E_up, but also progressively smaller portions of the envelope, due to the decreasing angular beam size (Tables 1 and 2). In cases such as the envelope of AFGL 2591, which is internally heated and enjoys a radially declining temperature profile, increasing J level and decreasing beam size both work in the same direction: the line signal probes progressively warmer gas. The observed velocity trend could be explained by (i) the very warmest core of the envelope being at a slightly different V_lsr than the surrounding envelope; or (ii) colder outer parts of the envelope being (mildly) influenced by blueshifted outflow material. The latter explanation is corroborated by the ¹³CO and C¹⁸O lines up to 3–2, which have redder V_centr than their 5 − 4 lines observed by HIFI (Fig. 1). The lines at E_up < 50 K lines are observed at IRAM 30 m or JCMT with beam sizes of 22′′ (¹³CO and C¹⁸O 1–0), 11′′ (C¹⁸O 2–1) and 15′′ (¹³CO 3–2 and HCO⁺ 5–4), all much smaller than the 38′′ beam applicable for the lowest-J HIFI lines.

To test the possible influence of the beam size on the centroid velocity value, we take advantage of the 1′  ×  1′ map of C¹⁸O 2–1 obtained at the IRAM 30 m telescope. When the map is convolved to simulate the largest Herschel beam, the velocity centroid is measured to move from −5.60 to −5.85 km s^-1, demonstrating that probing a larger portion of the envelope results in bluer line profiles. We conclude that the trend of increasing V_lsr with increasing J-level is likely to be caused by the systematically decreasing Herschel beam size (see Fig. 6).

Fig. 4 Spectral line energy distributions of CO a); its isotopic variants (b), c), d)); HCO+e); CS f); HCN g) and HNC h), as a function of Eup. Error bars on the observations (in black) represent the quadratic sum of the absolute intensity calibration uncertainty and the formal error on the Gaussian fit. Observed intensities represent only the envelope component for those lines fitted with a double Gaussian profile (Table 2). Since the intensity for a line fit by a single Gaussian is likely to be contaminated by outflow contribution as well, gray error bars provide an indication of the line intensity when up to half of the total line strength is attributed to the outflow. Downward pointing triangles indicate upper limits for undetected lines (see text in Sect. 3.1). Integrated line intensity values from the spherical model are marked by crosses and connected by dashed lines. Intensities from ground-based observations (J = 1–0, 2–1 and 3–2 for CO species, lowest transitions shown for HCO+, CS, HCN and HNC) and the corresponding model points (left of the dotted vertical line) are scaled down by a factor 10, as indicated by the “ ÷ 10” labels. Gray crosses accompanied by upward pointing arrows denote cases where the modeled line intensity falls off the vertical scale.

In addition, we investigate the trend of FWHM of the lines as a function of E_up in Fig. 5 (bottom), as measured in the envelope component of the Gaussian fits. Lines of ¹³CO and HCO⁺ are systematically narrower with increasing E_up and, again, decreasing beam size. As with the V_centr trend, the two ground-based lines for ¹³CO form an exception. The systematic shift is less clear or even reversed for C¹⁸O.

To explain systematically decreasing line widths, a turbulent core model (McKee & Tan 2003) could be invoked, in which turbulent velocity increases with radial distance due to energy dissipating from the warm inner regions into the colder outer envelope. For example, Herpin et al. (2012) successfully model an extensive set of H₂O lines with such a model description. If we adopt a similar model, however, line widths would be expected to increase toward lower transitions for all species. Since the observed FWHM values for C¹⁸O between 0 and 150 K and ¹³CO below 50 K show no such relation, we fall back on the hypothesis that the outflow material contributes to the “envelope component” line widths. As illustrated in Fig. 6, the smaller telescope beam at higher frequencies (and higher E_up) picks up less emission from the spatially extended outflow.

In addition, systematic trends of line velocity and line width with energy and/or beam size might be explained by an infall velocity gradient. This does not seem to be the case, however, for AFGL 2591, because, even at a spectral resolution of 0.1−1 km s^-1, none of our line shapes (Figs. 1, 2) show asymmetries that could indicate infall.

In conclusion, the trends in both V_lsr and FWHM suggest that the “envelope component” is kinematically affected by the outflow material, even after separating the obvious contribution of outflow material by means of the double Gaussian fits. This implies that our method to separate the outflow from the envelope is only partially successful. The assumption that the two physical components, quiescent envelope and outflow, are fully disconnected and do not influence each other dynamically obviously constitutes a highly simplified description of reality.

3.4. Angular size of emission region

Information about the angular size of an emission region can be derived from direct comparison of line intensities measured at the same frequency with different telescopes and beam sizes. One of the few atmospheric transmission windows that allows far-infrared frequencies to be observed from the ground is the region around 800 GHz, where this work presents the HCN 9–8 lines at 797.433 GHz and the CO 7–6 line at 806.652 GHz in a 27′′Herschel beam. These two lines were measured toward AFGL 2591 by Boonman et al. (2001) using the JCMT, with a beamsize of 8′′ (FWHM). The angular size of the emission region (θ_s) is coupled with the measured line intensities through its relation with intrinsic line brightness T_s, observed main beam line brightness T_mb, and telescope beam size θ_beam, (Wilson et al. 2009): (1)assuming that θ_s is small compared to the telescope beam. Because T_s should be independent of the observing facility, θ_s can be derived by substituting sets of (T_mb, θ_beam) for both JCMT and Herschel and equating the results.

We take the peak T_mb values of CO 7–6, HCN 9–8 and HCN (ν₂ = 1) 9–8 as measured at the JCMT, with an uncertainty of 50% (Boonman et al. 2001). In Table 3 we compare them to the peak intensities of the same lines measured in our Herschel/HIFI spectrum. For the vibrationally excited HCN, the HIFI line strength is an upper limit. Equation (1) is then used to derive limits on the angular size of the emitting region, θ_s. In the limiting case of point-like emission regions in both beams, i.e., θ_s² ≪ θ_beam², the line intensity ratio JCMT/Herschel would be 27²/8² ≈ 11, from simple beam dilution. For larger emission regions, the intensity ratio decreases.

The values for θ_s listed in Table 3 imply that the HCN 9–8 emission originates in a very confined region of the envelope of at most 1′′ in diameter. The CO 7–6 emission, however, is found to have an emission region of size θ_s ≳ 10′′. Although the latter violates the original assumption that θ_s ≪ θ_beam, it does indicate that CO 7–6 stems from a much more extended portion of the envelope than HCN 9–8. Implications of the sizes of emission regions derived here are discussed in Sect. 6.4.

Fig. 5 Trends of fitted centroid velocity (top) and fitted FWHM (bottom) of the envelope component versus Eup.

Fig. 6 Indication of the half-power beam width of Herschel/HIFI, represented by black circles: 42′′ at 500 GHz (left) and 19′′ at 1100 GHz (right). The contour map represents HCO+ 4–3 emission toward AFGL 2591 measured by the HARP-B instrument at JCMT. The blue- and redshifted wings are shown in solid and dashed contours, respectively, and are defined as in Van der Wiel et al. (2011). The lowest contour and spacing between consecutive contours is 5 K km s-1 for both components.

4. Envelope spectral line energy distribution model

4.1. Model setup

Radiative transfer calculations with the Monte Carlo ratran code (Hogerheijde & Van der Tak 2000) are employed to provide a model comparison to the observed SLEDs from Sect. 3.2. We rely on spectroscopic parameters and collisional rates from the lamda database (Schöier et al. 2005), which uses data provided by Yang et al. (2010) for CO, Dumouchel et al. (2010) for HCN and HNC, Flower (1999) for HCO⁺ and Turner et al. (1992) for CS.

The physical model we use is derived from a fit to the radial intensity profile of the sectors of the 450 and 850 μm JCMT SCUBA maps that appear spherical and to the λ ≳ 50 μm section of the continuum spectral energy distribution (Luis Chavarría, priv. comm.), with a strategy similar to that of Van der Tak et al. (1999, 2000). Our current model takes into account the new distance of 3.3 kpc and a corresponding input luminosity of 2.1  × 10⁵ L_⊙. Briefly, the model consists of a radial density gradient for molecular hydrogen described by an r^-1 power law (4  × 10⁷ down to 7  × 10² cm^-3), a total gas mass of ~500 M_⊙, and a corresponding temperature profile which varies monotonically from ~1400 K in the innermost shell to 27 K in the outer shell at 70 kAU. The radial distance at which the signal in the SCUBA 450 μm map falls below 3σ determines the value of the outer radius of the model envelope. The total luminosity, as determined by integrating over the best fit model continuum spectral energy distribution, is 8  × 10⁴ L_⊙; the difference with the observed luminosity is due to the absence of a cavity and a disk in our one-dimensional model, which does not include near-infrared emission. Jiménez-Serra et al. (2012) have found chemical segregation at very small scales (≲1′′, 3 kAU), but this is not expected to impact our models at the scales probed by our single dish observations (≳15′′, 50 kAU). Therefore, molecular abundances in our models are kept constant throughout the envelope and are scaled to match the observed SLEDs.

After the Monte Carlo method solves for the molecular level populations, ray tracing determines the emergent line intensity, which is convolved with a Gaussian telescope beam with the appropriate HPBW at the relevant frequency. Resulting spectral maps are collapsed along the spectral direction to obtain the integrated line intensity value, which is compared to observations in Fig. 4. Because the outflow component is not represented in the model, we compare predicted integrated line intensities from the models described above only with the envelope component of the observed spectral lines (Sect. 3.2).

Despite the obvious increase in linear size scales when comparing an older model for AFGL 2591 assuming a distance of 1 kpc (e.g., Van der Tak et al. 2000; Van der Wiel et al. 2011) with the one employed in this work, the differences in emergent line intensities are relatively minor. This is due to the similar temperature and density structure, which lead to comparable excitation conditions. A noteworthy element of the new model is the temperature of >1000 K in the central shells of the envelope, compared to 372 K for the old model. The maximum temperature in the model is essentially set by the arbitrary choices of the thickness of the spherical model shells and the size of a central sphere which is vacant of gas. By running test models in which the very hottest inner shells are removed, we determine that the hot gas in the central envelope has no impact whatsoever on the emergent line intensities of rotational levels up to E_up = 700 K. Essentially, none of the line emission emanating from the hot gas escapes the envelope to be detected by external observers.

4.2. Comparison of model results with observations

With a CO abundance of 5  × 10^-5 relative to H₂, one quarter of the 2  × 10^-4 applicable for warm gas (Van der Tak et al. 1999), the modeled SLED of the main form of CO lies a factor of 2 to 4 above the observed points up to J = 9–8 (off the scale in Fig. 4a). Beyond that, the model SLED drops off sharply and underproduces observed line intensities from J = 13–12 onward by factors of 2 to 3. Although with smaller discrepancies, the models for the isotopic variants ¹³CO, C¹⁸O and C¹⁷O (with abundances of 5  × 10^-7, 6  × 10^-8 and 3  × 10^-8, respectively) show a similar trend in Figs. 4b–d: low-J emission is over-predicted while high-J emission is under-predicted.

The HCO⁺ model, at an abundance of 1  × 10^-8, seems successful: it matches the envelope component of the observed lines J = 6–5 and higher to within a factor 2 (Fig. 4e). However, the same model overproduces HCO⁺ 4–3 by a factor 5. The models for the other species studied here, CS, HCN and HNC (abundances of 4, 6 and 1  ×  10^-8, respectively), again overproduce low-J emission by factors ~5 and are marginally consistent with or fall below the observed points at higher J (Figs. 4f–h).

The general picture that emerges is that our smooth, spherically symmetric models predict too much emission from cold gas (~50 K) or too little emission from warm gas (≳100 K). This effect persists even when we consider the possibility that the envelope contribution to the emission is overestimated by a factor of 2 for lines that are fitted by a single Gaussian (gray error bars in Fig. 4). For all species studied here, choosing a different molecular abundance has an effect on the vertical offset of the model SLED, especially at more optically thin transitions at high energies, but it does not sufficiently reduce the slope between E_up ~ 100 and 200 K to match the observations. Our model serves to show that a spherical, homogeneous geometry offers a poor match to the true excitation balance. Therefore, the absolute molecular abundances used in this section are not meaningful by themselves. For example, the CO abundance in our model is lower than generally assumed, but still severely overproduces low-J line intensities and optical depths.

A model might be considered where CO is depleted from the gas phase, to suppress specifically the low-J CO emission. However, the temperature profile we use for the envelope does not reach values below 25 K. Hence, the model gas is too warm for CO freeze-out to occur, consistent with the absence of CO ice features in mid-infrared spectra (Van der Tak et al. 1999).

Our model is based on a density profile with a power law index of α = −1.0, as outlined in Sect. 4.1 and also used in earlier work on this source (see Sect. 1). We have explored ratran models with steeper profiles (α = −1.5), but the resulting SLEDs are very similar in shape to the ones with α = −1.0. A steeper power law therefore does not mitigate the mismatch of the shape of the linear rotor SLEDs with the observations.

Considering the spectral dimension, the modeled and beam-convolved line profiles of, e.g., CO and HCO⁺ up to 7–6 and HCN up to J = 10–9 appear broadened or even self-absorbed, indicating optically thick emission. For example, the optical depths at line center determined during the ray tracing, for the main isotopic form of CO, are between 90 and 10 for 3–2 up to 11–10, while only the 16–15 is marginally optically thin at 0.8. The modeled lines for the isotopologues of CO are all optically thin, as are all CS transitions except J = 7–6. However, HCO⁺ is optically thick until 10–9, HCN until 9–8 and HNC until 7–6. Indications for such high optical depths are not seen in the observed line profiles (Figs. 1 and 2). This discrepancy between observed and modeled line shapes confirms earlier suspicions (Van der Tak et al. 1999; Van der Wiel et al. 2011) that line optical depths are too high in a spherical model. In Sect. 6.1, we will suggest possible approaches to address the optical depth with more sophisticated models.

Fig. 7 Zoom of the observed 0 km s-1 feature in various CO lines (histograms). The results from the radex slab are shown as solid lines. To improve visibility, the brightnesses for C18O (both observations and models) are multiplied by 10. Vertical offsets are applied as follows: +5 K for CO 5–4, −5 K for C18O 2–1, −10 K for C18O 1–0 and −15 K for 13CO 1–0. The gray dashed line marks the approximate centroid velocity of the absorption feature at + 0.2 km s-1. The model representations only include a component for the main envelope emission for 5–4 and 3–2, where a background is needed to absorb against. The 2–1 and 1–0 model lines only account for the foreground emission at ~0 km s-1, hence the discrepancies at Vlsr ≲ −1 km s-1 where the main envelope of AFGL 2591 dominates.

5. Physical conditions per component

5.1. Foreground clouds

Our CO 5–4 spectrum exhibits an absorption feature near 0 km s^-1 (Fig. 7), which is also seen in JCMT data of CO 3–2 and is coincident in velocity space with a known foreground cloud (Van der Tak et al. 1999). Conversely, an emission peak is visible near 0 km s^-1 in our IRAM spectra of ¹³CO 1–0 and C¹⁸O 1–0 (E_up = 5.3 K, Fig. 7). Mitchell et al. (1992) show spectra of CO 2–1 toward various positions in the envelope showing the same absorption feature near V_lsr = 0 km s^-1. It is not seen in higher-J lines of CO, indicating that the absorbing foreground cloud must be relatively cold. It must also be of relatively low column density, given that the absorption feature is not present in ¹³CO and C¹⁸O 5–4.

Minh & Yang (2008) observe an extended 0 km s^-1 component in CO and ¹³CO 1–0 emission, peaking 8′ to the northeast and northwest of AFGL 2591. Our observations show that this gas cloud also extends to the position of AFGL 2591 itself, and that it must be in the foreground, since it absorbs CO 3–2 and 5–4 emission from the warm protostellar envelope.

We construct a simple “slab” component to account for the foreground CO absorption, using the interface to radex (Van der Tak et al. 2007) available in cassis. First, the shape of the CO 5–4 emission line is mimicked by three “LTE” components with excitation temperatures of 100 K, V_centr of −4.5, −7.3 and +1.6 km s^-1, FWHM values of 5.1, 15.7 and 2.5 km  s^-1 and CO column densities of 6, 9 and 0.3  × 10¹⁸ cm^-2. Note that these emission components are merely invoked to create representative background emission and are not meant to represent real physical conditions in the protostellar envelope. On top of the mock background emission, a radex slab is added with N_CO = 3 × 10¹⁷ cm^-2 and a kinetic gas temperature (T_kin) of 10 K. This choice for temperature is motivated by the fact that it should be lower than the minimum gas temperatures of ~25 K reached at the outer edges of the protostellar envelope of AFGL 2591, while at the same time it should be high enough to result in detectable emission in the 1–0 transition of ¹³CO and C¹⁸O, with E_up = 5.3 K.

Because the foreground component is likely to be part of a larger scale cloud (Minh & Yang 2008), the angular size of the slab is set at 50′′, larger than any of the beam sizes of our observations. Finally, the absorption component is given an intrinsic FWHM of 0.8 km s^-1 and V_lsr of +0.3 km s^-1 to match the width and position of the observed absorption feature. After the frequency-dependent optical depth calculation and the radiative transfer, the result is convolved with the relevant telescope beam to make a direct comparison to observed line profiles in Fig. 7.

With the parameters for the cold material chosen above to match CO 5–4, we also find self-absorption in CO 3–2 (again, against a mock background emission line) that is qualitatively consistent with the JCMT observations. Moreover, due to the 10 K gas temperature, the same component naturally explains the emission bumps observed in ¹³CO 1–0, C¹⁸O 1–0 and even in C¹⁸O 2–1 (Fig. 7). The column densities needed for ¹³CO and C¹⁸O are 3  ×  10¹⁵ and 2 × 10¹⁴ cm^-2, respectively.

We conclude that a cold slab of foreground gas explains the 0 km s^-1 feature. Assuming an abundance of CO/H₂ of at most 10^-4, we derive a column density N_H2 ≳ 3 × 10²¹ cm^-2. The visual extinction (A_V) for this column density would be ~3 mag, following Bohlin et al. (1978), purely for the molecular part. If we add to this 1 mag each (front and back) for the atomic surfaces of the cloud where CO is photodissociated, we find that the total A_V ≳ 5. The derived conditions in the 0 km s^-1 component (“foreground A”) are summarized in Table 4.

Whereas the 1–0 emission features are insensitive to changes in T_kin between 5 and 20 K, the depth and shape of the CO 5–4 absorption features constrains T_kin to [9,15] K. The value adopted for N_CO is motivated by the depth of the self-absorption feature, which starts to deviate significantly when N_CO is changed by a factor ≳2. The model result is insensitive to the number density of the collision partner in radex as long as n(H₂) is above ~10⁴ cm^-3, i.e., when the low-J CO excitation balance is in LTE.

Given the inferred nature of the foreground cloud, it is important to consider that H₂ densities in translucent clouds are generally in the order of 10³ rather than >10⁴ cm^-3. With a lower H₂ density, the J = 1–0 emission features in ¹³CO and C¹⁸O can be explained with temperatures as high as 17 K. Hence, the density as well as the upper limit to the kinetic temperature of the foreground cloud are poorly constrained.

In addition, a discrepancy is found with respect to the observed CO 6–5 line profile that shows no absorption at 0 km s^-1. In fact, if any feature is present at that velocity, it is in emission (Fig. 1). The slab model with T_kin = 10 K and N_CO = 3 × 10¹⁷ cm^-2, however, predicts a detectable absorption feature for CO 6–5. We are unable to explain the CO 6–5 without compromising the match to the lower-J 0 km s^-1 feature.

Although not included in the considerations in this paper, several additional probes of the 0 km s^-1 foreground cloud are available. Firstly, besides the CO signatures shown in Fig. 7, the ³P₁–³P₀ line of neutral carbon (at 492.161 GHz in our CHESS spectral scan) also shows a weak emission feature near 0 km s^-1. Secondly, an upper limit on the density might be derived from the fact that the 0 km s^-1 component is not seen in lines of, e.g., HCO⁺ and HCN (Fig. 2). Thirdly, signatures of CH, CH⁺, H₂O⁺ and other hydride species do show absorption near 0 km s^-1 (Bruderer et al. 2010b), as does the HF signature in the CHESS spectrum detected by Emprechtinger et al. (2012). A more extensive exploration incorporating all observed tracers of the 0 km s^-1 component toward AFGL 2591 is therefore warranted.

In addition to the 0 km s^-1 component described above, there is another absorption component, at +13 km s^-1, visible in Herschel observations of HF (Emprechtinger et al. 2012, but not discussed there) as well as H₂O (Choi et al., in prep.). It is not seen, however, in any lines of CO species presented in this work, neither those observed with HIFI (Fig. 1) nor those observed at IRAM 30 m and JCMT (the velocity range covered in the observations is larger than what is shown in Fig. 7). The diatomic hydrides observed in the Herschel WISH program (Bruderer et al. 2010b) do not show clear evidence of emission or absorption at + 13 km s^-1 either.

Based on our ¹³CO 1–0 spectrum, the + 13 km s^-1 feature would need to have a peak intensity of ≲ 0.1 K to ensure a non-detection. With that upper limit, a brief exploration of radex models tells us that the ¹³CO column density is at most 10¹⁴ cm^-2. Using the canonical scaling factors to convert to CO (60) and noting that the CO abundance relative to H₂ can be as low as 10^-6 in diffuse gas, we derive that N_H2 ≲ 6  × 10²¹ cm^-2 in the + 13 km s^-1 component. In addition, the CHESS spectrum of HF 1–0 presented by Emprechtinger et al. (2012) provides a lower limit for the column density of the + 13 km s^-1 cloud. While Emprechtinger et al. do not discuss the + 13 km s^-1 component, an absorption feature does appear at that velocity with a depth of ~70% relative to the saturated 0 km s^-1 component. For the latter, they derive N_HF > 4  × 10¹³ cm^-2. Scaling by the relative depth of the absorption features and assuming that all fluorine is in the form of HF (HF/H₂ = 3.6  × 10^-8), we find that N_H2 ≳ 8  × 10²⁰ cm^-2. Column density limits are summarized in Table 4 (“foreground B”).

Fig. 8 Lines of constant intensity ratio of C18O 5–4/7–6 (solid lines) and HCO+ 4–3/7–6 (dashed lines), resulting from a radex grid calculation. Curves are labeled by their line ratio value. The gray areas mark the 1-σ uncertainty ranges of the observed line intensity ratios in the outflow component, including the fit error (see Table 2) and the 10% calibration uncertainty (see Sect. 2.1). The darker gray area denotes the region where the observed line intensity ratios of both tracers overlap.

5.2. Outflow gas

Previous studies indicate that CO species, HCO⁺, and neutral atomic carbon are present in the outflow in significant amounts (e.g., Choi et al. 1994; Hasegawa & Mitchell 1995; Doty et al. 2002). Our results confirm that line profiles of CO species and of HCO⁺ are well represented by two separate kinematic components (Sect. 3.2).

Using radiative transfer calculations with the radex code, we estimate kinetic temperatures, number densities and column densities of CO and HCO⁺ in the outflow component. This approach requires an assumption of an isothermal homogeneous material, in which the radiative transfer is solved using collisional rates from the LAMDA database (Schöier et al. 2005). Following Van der Tak et al. (2007), we construct grids of radex models and calculate line intensity ratios for values of n_H2 between 10³ and 10⁸ cm^-3 and T_kin between 10 and 400 K. In the input for the grid calculation, the column densities are intentionally low in order to keep all modeled lines optically thin. In this limit, the intensity of each line is directly proportional to the column density and line ratios are therefore insensitive to the absolute column density. In Fig. 8 we show the resulting line ratios of HCO⁺ 4–3/7–6 and C¹⁸O 5–4/7–6 as a function of n_H2 and T_kin. When the observed line intensity ratio is indicated in the diagram, it is possible to constrain a region in parameter space that fits the observations.

First, the mere presence of high-J HCO⁺ emission in the outflow component is an indication that number densities must be relatively high, given that critical densities are in the order of 10⁷ cm^-3 for the 6–5, 7–6 and 8–7 transitions. Second, from the C¹⁸O line ratio in Fig. 8, it is evident that densities below 10⁴ cm^-3 require temperatures that are unphysically high. The observed HCO⁺ and C¹⁸O line ratios intersect in a region in (n_H2,T_kin)-space near 10^4.7–10⁶ cm^-3 and ~60–200 K. Inside this region, higher densities imply lower temperatures, and vice versa (Fig. 8). The measured line ratio for ¹³CO 5–4/7–6 gives the same constraints as C¹⁸O, confirming that signal-to-noise and optical depth do not affect our result. Finally, a line ratio diagram for C¹⁸O 1–0/2–1 suggests the presence of an additional component of colder gas at T_kin < 50 K. However, this result should be considered with caution, because (i) the kinematic separation of envelope and outflow contributions is not unambiguous (Sect. 3.2); and (ii) the beam filling factor of the outflow component could be different for the 2–1 and 1–0 observations (beam HPBW 11′′ and 22′′, respectively).

With n_H2 = [10⁵,10⁶] cm^-2 and T_kin = [60, 200] K, a column density in the order of 10¹³ to 10¹⁴ cm^-2 for HCO⁺ reproduces the line strengths observed for the outflow. Similarly, a C¹⁸O column density of 3 × 10¹⁵ cm^-2 ( ± 30%) is required to explain the HIFI C¹⁸O line strengths in the outflow component. With a CO/C¹⁸O ratio of 500 and a CO abundance of 10^-4 relative to H₂, we derive N_H2 = 1.5 × 10²² cm^-2 and a molecular abundance in the order of 10^-9 to 10^-8 for HCO⁺ in the outflow gas. This HCO⁺ abundance falls in the same range as previously adopted values for the envelope, ranging from 4 × 10^-10 (Carr et al. 1995) to ~1  ×  10^-8 (Van der Wiel et al. 2011).

In conclusion, Table 4 contrasts the physical conditions derived here for the protostellar outflow component (Sect. 5.2) with those found for the CO foreground cloud (Sect. 5.1) and in the envelope. Our current analysis indicates that the outflow gas is not significantly different from that in the envelope, considering gas density, gas temperature, as well as the chemical balance of CO and HCO⁺.

6. Discussion and conclusions

6.1. Envelope structure

The spherically symmetric radiative transfer models are unable to simultaneously explain the low- and high-J transitions of the observed integrated line intensities (Sect. 4, Fig. 4): the modeled SLEDs of all species studied here decline too steeply between ~100 and 200 K. This indicates that the passively heated, spherical envelope structure does not populate the high J levels sufficiently, or that emission from those levels is too optically thick. Line optical depths are known to be higher than observed for AFGL 2591 in ground-based low excitation lines (Van der Wiel et al. 2011) and in Sect. 4 we have noted that optical depths for the modeled higher J lines are lower, but still >1 in many cases.

To boost line intensities at E_up above ~150 K without increasing the lower level intensities too much, it would be appropriate to apply a model geometry that includes outflow cavities. This would affect not only the density structure to create lower opacity pathways (if a suitable inclination angle is chosen), but also the chemical and thermal balance through direct illumination of envelope material by UV radiation from the central massive protostar. Indeed, we demonstrate in Sect. 4 that the hot gas (>1000 K) in the center of the model envelope does not contribute to emergent line intensity if it is embedded in an isotropic envelope.

In this light, our AFGL 2591 observations should be compared with a model scheme such as that by Bruderer et al. (2009a, 2010a, 2012), with a physical geometry appropriate for the massive protostellar envelope of AFGL 2591. Their model allows for an outflow cavity to be present in an axisymmetric (or even three-dimensional) density distribution, from which the dust temperature profile is derived through continuum radiative transfer. The UV field is calculated and used as input to iteratively determine the position-dependent chemical balance, excitation balance, gas temperature and resulting line emission. The fully self-consistent chemical treatment includes freeze-out, photodissociation, X-ray ionization and other processes. Preliminary exploration of similar models already shows that, as the cavity allows a larger volume of lower density gas to be irradiated and heated directly, SLED shapes become significantly flatter compared with spherically symmetric geometries (Van der Wiel 2011, Chapter 4; Steven Doty 2011, priv. comm.). Bruderer et al. (2012) also show that in exposed layers of the gas envelope, the gas and dust temperature become decoupled, rendering inadequate our model for which the gas temperature is derived from dust continuum measurements (Sect. 4).

The HCO⁺ molecule in particular could provide an interesting probe of outflow cavity geometry, since, contrary to CO and HCN, the HCO⁺ abundance is expected to be enhanced in the cavity walls (Simon Bruderer 2011, priv. comm.). Additional constraints could come from the incorporation of CO⁺ observations. Although no signatures from this species are detected in the spectral scans with Herschel/HIFI (Paper II) and with JCMT/HARP-B (Van der Wiel 2011, Sect. 3.A), two submillimeter CO⁺ lines are detected in deeper spectra by Stäuber et al. (2007). From the sample of low-mass and high-mass protostars in their study, they conclude that X-rays play a role in the dissociation and excitation of molecules in low-mass protostellar envelopes, but that CO⁺ observations can be explained by far-UV radiation in high-mass objects. Even in dissociation regions, CO⁺ emission could be weak or absent, since this ion is easily destroyed by collisions (Stäuber & Bruderer 2009).

Additionally, line opacity properties of the modeled protostellar envelope can be adjusted by distributing the same amount of gaseous material inhomogeneously, as previously suggested by Van der Wiel et al. (2011). This could be used to further alleviate the optical depth effects (Sect. 4), if the introduction of cavities in a two-dimensional geometry does not prove sufficiently effective to match the observations.

Studies of low-mass star-forming envelopes have shown that a passively heated envelope is the dominant component for lines with E_up ≲ 200 K, but that the ~200–4000 K section of the CO SLED is a discriminating probe of the contribution of UV irradiation and shocks (Van Kempen et al. 2010; Visser et al. 2012). While the highest J CO line presented in this work is 16–15, the highest attainable to HIFI, Herschel/PACS has observed AFGL 2591 at frequencies covering lines up to J_up ≈ 30, to be analyzed in the context of the WISH key program (Van Dishoeck et al. 2011). We expect the extension of the CO ladder to higher energies to shed more light on excitation drivers in massive protostellar envelopes such as AFGL 2591.

6.2. Foreground clouds

We identify a 0 km s^-1 foreground cloud in Sect. 5.1, for which we derive a kinetic gas temperature between 9 and ~17 K, based on the combination of absorption and emission signatures in the J = 1–0 up to 5–4 transitions of CO and isotopologues. The material in the foreground cloud discussed here is redshifted with respect to AFGL 2591. Therefore, the foreground could be feeding additional gas to the protostellar system, but evidence would have to be found to support physical proximity of the 0 km s^-1 cloud to AFGL 2591.

For a second foreground cloud, at V_lsr = +13 km s^-1, we use our non-detections in CO and the absorption in HF to derive limits for the column density: 8  × 10²⁰ cm^-2 ≲ N_H2 ≲ 6 × 10²¹ cm^-2 (Sect. 5.1). The CO abundance relative to H₂ must be ≲10^-5 for the derived limits to be consistent with each other (for CO/H₂ = 10^-4, the upper limit would be 6  ×  10¹⁹ cm^-2). In such diffuse conditions HF is more abundant than CO, as predicted by Neufeld et al. (2005) and observed, for example, by Sonnentrucker et al. (2010) on the line of sight to W51. It is surprising, however, to detect a relatively strong absorption signature of H₂O without the presence of CO. Water abundance values will be derived by Choi et al. (in prep.).

6.3. Outflow gas

The kinetic temperature of 60–200 K and number density of ~10⁵–10⁶ cm^-3 derived for the outflow gas are comparable to those in the passive envelope (Sect. 5.2). The component that we identify a few km s^-1 blueward of the systemic V_lsr is kinematically distinct from the envelope (Table 2), but it apparently harbors gas that is similar in density and temperature. Our results are consistent with conditions found in other protostellar outflows (e.g., Giannini et al. 2001; Moro-Martín et al. 2001; Lefloch et al. 2010; Yıldız et al. 2010, 2012; Bjerkeli et al. 2011; Gómez-Ruiz et al. 2012). Specifically, CHESS spectra toward the L1157-B1 outflow reveal H₂ gas densities of ~10⁵ to a few 10⁶ cm^-3 and temperatures around 200 K (Codella et al. 2012; Benedettini et al. 2012; Lefloch et al. 2010, 2012).

We speculate that the similar physical and chemical conditions (Table 4) between outflow and envelope could point to outflow gas that has been expelled from the envelope relatively recently. The blue component observed in this work may in fact trace swept-up envelope material ( ≲ 10 km s^-1 from the systemic velocity) rather than fast outflowing gas that is farther separated from the parent envelope, see, e.g., Van der Tak et al. (1999) for AFGL 2591, or Qiu et al. (2011) and Gómez-Ruiz et al. (2012) for other protostars. Our outflow conditions are compatible with the temperature and density derived for the warmest of the three outflow components in L1157-B1, as found by Lefloch et al. (2012). These authors argue that lower kinetic temperatures indicate older outflows, analogous to our conclusion that the warm, dense outflow material from AFGL 2591 is young. We should also caution that the double Gaussian fits that are applied to distinguish the components (Sect. 3.2) are likely to suffer from cross-contamination. Both components could partly trace the same gas and, as a consequence, the derived physical conditions for the envelope and outflow may be artificially equalized.

6.4. Compact HCN emission

In Sect. 3.4, our observations of HCN 9–8 and CO 7–6 are compared with those from a larger telescope, in order to impose limits on the size of the emitting regions. We find sizes of ≲ 1″ for HCN 9–8 and ≳ 10′′ for CO 7–6. Benz et al. (2007) find source sizes of < 1′′ in an interferometric HCN 4–3 map, but the extended emission in their observation is likely filtered out. Conversely, single dish observations by Van der Wiel et al. (2011) reveal larger scale (>30″) emission in the 4 − 3 transition of HCN, but this line traces lower temperatures and lower densities than HCN 9–8. Considering our envelope model from Sect. 4, only in the very inner shells (<30 AU ≈ 0.01′′) does the number density come within an order of magnitude of the critical density of HCN 9–8, 3  × 10⁸ cm^-3. It is therefore conceivable, and fully consistent with our observations, that the HCN J = 9–8 emission arises in a region smaller than 1′′ in angular size. The critical density for CO 7–6 is a factor ~50 lower than that of HCN 9–8; it is therefore expected that CO 7–6 emerges from a more extended region of the envelope.

7. Summary and outlook

In this first Herschel CHESS paper on the massive star-forming envelope AFGL 2591, we present HIFI spectral line measurements of a wide range of transitions (E_up/k = 80 up to 300 − 700 K) of CO, ¹³CO, C¹⁸O, C¹⁷O, HCO⁺, CS, HCN and HNC, supported by ground-based measurement of lower excitation levels. Emission lines are separated into a narrow “envelope” and a broader “outflow” component, if the S/N is sufficient (Sect. 3.2). The integrated line intensities of the envelope component, presented in the form of spectral line energy distributions (SLEDs), are compared with radiative transfer results based on a spherically symmetric envelope model (Sect. 4).

The main conclusions of this work are the following.

• 1.

  The model SLEDs for all species fall off too steeply between100 and 200 K to explain the flat shape of theobserved envelope SLEDs (Sect. 4). We expectthat a model including outflow cavities and direct UVillumination will provide a better match(Sect. 6.1).

• 2.

  The outflow gas is found to be as dense and warm as that in the quiescent envelope (Sect. 5.2). We conclude that the outflow gas traced here may have been expelled only recently (Sect. 6.3).

• 3.

  Trends of line centroid and line width with frequency, excitation level or telescope beam size are identified in the envelope component. In absence of a single, consistent, physical explanation, we ascribe these trends to ambiguities in the separation of the envelope and outflow components (Sect. 3.3).

• 4.

  By comparing line intensities measured with Herschel (3.5 m) and JCMT (15 m), we conclude that HCN 9–8 emission must originate from a region < 1′′ in angular size (Sects. 3.4 and 6.4).

• 5.

  A known foreground cloud at V_lsr = 0 km s^-1 is detected in several lines of CO and isotopologues. Limits to its temperature and column density are derived in Sect. 5.1.

• 6.

  A second foreground cloud, detected in HF and H₂O at + 13 km s^-1, remains undetected in our observations (Sect. 5.1). Its nature and column density are discussed in Sect. 6.2.

Future papers in this series based on the CHESS observations of AFGL 2591 will be ordered by chemical families of species, e.g., N-bearing species (Paper II), S-bearing species. Paper II will also include details of all other molecular line signal detected in the 490–1240 GHz spectral scan. In addition, future work should be devoted to the application of a more sophisticated model geometry for AFGL 2591 to explain the highly excited rotational lines accessible with Herschel. Such a model should include at the least an outflow cavity geometry and chemical effects of direct UV illumination (now available for hydrides and C/C⁺/CO, Bruderer et al. 2010a,b, 2012), but possibly also shock physics and inhomogeneous envelope structure.

Acknowledgments

The authors acknowledge constructive discussions with members of the HIFI Instrument Control Center. We thank Mihkel Kama and the other members of the CHESS data reduction team for frequent interaction on HIFI spectral scan data processing. We thank Luis Chavarría for providing the new spherical envelope model structure, and Marco Spaans, Steve Doty, René Plume and Simon Bruderer for discussions on physical and chemical modeling. We acknowledge constructive comments by an anonymous referee and by the Editor. 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 US. 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 Zürich, FHNW; USA: Caltech, JPL, NHSC. This work is supported in part by the Canadian Space Agency and NSERC.

References

All Tables

All Figures

Fig. 1 Spectral line profiles of CO, 13CO, C18O, and C17O. The 1–0 and 2–1 spectra are from the IRAM 30 m telescope, the 3–2 lines are from JCMT/HARP-B. All higher-J lines are from Herschel/HIFI. For clarity, consecutive spectra, each labeled by their upper and lower J-level, are offset in the vertical direction. Corresponding rest frequencies and energy levels are listed in Table 2. The systemic velocity ( − 5.8 km s-1) of AFGL 2591 is indicated by the dashed vertical line.

In the text

	Fig. 2 As Fig. 1, but for HCO+, CS, HCN, and HNC. All lines are measured with Herschel/HIFI, except HCO+, HCN, and HNC 4–3, as well as CS 7–6, which are measured with JCMT.
In the text

	Fig. 3 Examples of Gaussian component fits to HCN 6–5 (top, single component) and HCO+ 6–5 (bottom, two components). Observed line profiles are shown as histograms, fit components are shown as smooth solid curves and the sum of the two components in the bottom panel is represented by the dashed curve.
In the text

Fig. 4 Spectral line energy distributions of CO a); its isotopic variants (b), c), d)); HCO+e); CS f); HCN g) and HNC h), as a function of Eup. Error bars on the observations (in black) represent the quadratic sum of the absolute intensity calibration uncertainty and the formal error on the Gaussian fit. Observed intensities represent only the envelope component for those lines fitted with a double Gaussian profile (Table 2). Since the intensity for a line fit by a single Gaussian is likely to be contaminated by outflow contribution as well, gray error bars provide an indication of the line intensity when up to half of the total line strength is attributed to the outflow. Downward pointing triangles indicate upper limits for undetected lines (see text in Sect. 3.1). Integrated line intensity values from the spherical model are marked by crosses and connected by dashed lines. Intensities from ground-based observations (J = 1–0, 2–1 and 3–2 for CO species, lowest transitions shown for HCO+, CS, HCN and HNC) and the corresponding model points (left of the dotted vertical line) are scaled down by a factor 10, as indicated by the “ ÷ 10” labels. Gray crosses accompanied by upward pointing arrows denote cases where the modeled line intensity falls off the vertical scale.

In the text

	Fig. 5 Trends of fitted centroid velocity (top) and fitted FWHM (bottom) of the envelope component versus Eup.
In the text

Fig. 6 Indication of the half-power beam width of Herschel/HIFI, represented by black circles: 42′′ at 500 GHz (left) and 19′′ at 1100 GHz (right). The contour map represents HCO+ 4–3 emission toward AFGL 2591 measured by the HARP-B instrument at JCMT. The blue- and redshifted wings are shown in solid and dashed contours, respectively, and are defined as in Van der Wiel et al. (2011). The lowest contour and spacing between consecutive contours is 5 K km s-1 for both components.

In the text

Fig. 7 Zoom of the observed 0 km s-1 feature in various CO lines (histograms). The results from the radex slab are shown as solid lines. To improve visibility, the brightnesses for C18O (both observations and models) are multiplied by 10. Vertical offsets are applied as follows: +5 K for CO 5–4, −5 K for C18O 2–1, −10 K for C18O 1–0 and −15 K for 13CO 1–0. The gray dashed line marks the approximate centroid velocity of the absorption feature at + 0.2 km s-1. The model representations only include a component for the main envelope emission for 5–4 and 3–2, where a background is needed to absorb against. The 2–1 and 1–0 model lines only account for the foreground emission at ~0 km s-1, hence the discrepancies at Vlsr ≲ −1 km s-1 where the main envelope of AFGL 2591 dominates.

In the text

Fig. 8 Lines of constant intensity ratio of C18O 5–4/7–6 (solid lines) and HCO+ 4–3/7–6 (dashed lines), resulting from a radex grid calculation. Curves are labeled by their line ratio value. The gray areas mark the 1-σ uncertainty ranges of the observed line intensity ratios in the outflow component, including the fit error (see Table 2) and the 10% calibration uncertainty (see Sect. 2.1). The darker gray area denotes the region where the observed line intensity ratios of both tracers overlap.

In the text