© ESO, 2015

1. Introduction

Stars are born inside dense molecular structures in the interstellar medium (ISM) that consist of gas and dust. The dust grains in clouds associated with star-forming regions absorb the short–wavelength radiation from the central stars, heat up, and then re-emit radiation at far-IR and sub–millimeter wavelengths. The thermal radiation from dust is optically thin at these wavelengths and thus is a good tracer of such physical parameters as temperature, density, and the gas mass of the clouds. Dust is efficiently heated through near IR radiation thanks to the broadband absorption properties of the dust grains. On the other hand, the gas temperature T_gas is mainly governed by heating processes driven by UV radiation (photoelectric heating, H₂ excitation, dissociation).

In the well-shielded centers of massive cores, the primary effects for the thermal balance are cosmic ray heating, molecular line cooling, and collisional heating or cooling due to dust-gas collisions (Goldsmith & Langer 1978; Tielens 2005; Draine 2011).

At densities above ~10^4.5, the dust and gas are expected to be collisionally coupled and are characterized by the same temperature (Goldsmith 2001). At lower densities, the rate of collisions between dust and gas decreases, and the cooling of the gas via fine-structure and molecular rotational transitions becomes dominant (e.g. CO, [C II] , and [OI]; Sternberg & Dalgarno 1995; Kaufman et al. 1999; Meijerink & Spaans 2005).

Excitation of the rotational levels of molecules like CO (low dipole moment) observed at (sub)mm wavelengths is provided by collisions with H₂. If these collisions are frequent enough to exceed the spontaneous decay rate, the levels will get into thermal equilibrium with H₂. The critical densities n_cr of CO and isotopologues for the J = 1 − 0 and 2−1 transitions at 50 K are ~ 2 × 10³ cm^-3 and ~ 2 × 10⁴ cm^-3, respectively¹. At this point their excitation temperature will approach the gas kinetic temperature, so that the study of these lines is ideal for gas kinetic temperature estimates. The rarer isotopologues of CO, such as ¹³CO and C¹⁸O, can be used to probe regions of high column density where the lines of the abundant isotopologues become optically thick. The critical densities n_cr of CO and isotopologues at 50 K are ~1.2 × 10⁶ cm^-3 for the J = 9 − 8 and ~ 1.5 × 10⁶ cm^-3 for the 10−9 transition, respectively.

S 140 is a well studied HII region that lies at the southwest edge of the molecular cloud L 1204 at a distance of 746 pc (Hirota et al. 2008). A larger distance of 910 pc has been adopted in previous studies (e.g., Crampton & Fisher 1974; Dedes et al. 2010; Ikeda & Kitamura 2011). This edge is illuminated by the B0V star HD 211880, creating a visible HII region and a photon-dominated region (PDR) at the ionization front (IF). The cloud also hosts a cluster of embedded high-mass young stellar objects (YSOs), IRS 1 to 3 (Evans et al. 1989) at a projected distance of ~ 75′′ northeast of the IF. A large number of low-mass stars are also forming, making S 140 an ideal laboratory for studying both the star formation of various masses and the differences caused by irradiation from external (B0V star) and internal sources (IRS 1 to 3).

The analysis of physical models around IRS 1 has been performed by a number of previous authors including Harvey et al. (1978), Gürtler et al. (1991), Minchin et al. (1993), van der Tak et al. (2000), Mueller et al. (2002), de Wit et al. (2009), and Maud & Hoare (2013). Poelman & Spaans (2006) report the clumpiness towards S 140 and give a value of n_H2 ~ 10⁴ cm^-3 for the interclump medium and ~4 × 10⁵ for the clump gas. Previous models lacked spatial information at many wavelengths, and many of them were based on the assumption that the gas temperature is coupled to the dust as a result of the high densities (e.g., Poelman & Spaans 2006).

In this study we have analyzed the temperature and column density around the S 140 high-mass star–forming region of embedded young stellar objects by combining Herschel-PACS and HIFI data with ground-based mapping observations of mm–wave lines using the IRAM 30 m telescope. The new PACS continuum data and IRAM spectroscopic maps provide us with the highest angular resolution spatial information available, covering an area that includes the YSOs in both datasets. The IF was covered by our larger IRAM maps (4′× 4′) but not by the PACS observations, since the latter cover a smaller area (45′′× 45′′). Furthermore, the data give us the opportunity to study independently of the gas and dust temperatures around FIRS 1−3 and check whether the assumption of well-coupled dust and gas is valid throughout this cloud.

In the following sections, we present the dust and gas modeling, along with the resulting parameters. We present the observations in Sect. 2, the observational results in Sect. 3, the gas and dust analysis with a comparison between gas and dust in Sect. 4, the effect of density gradients using more advanced dust and gas models in Sect. 5, and a discussion of our results and conclusions in Sect. 6.

2. Observations and data reduction

2.1. PACS data

To investigate the temperature and column-density structure of the dust we analyzed a single footprint of PACS (Poglitsch et al. 2010) with its 5 × 5 spatial array of ~ 9′′ pixels over a range of wavelengths from ~ 70 μm to ~ 200 μm (obsId: 1342222256). The data was reduced using HIPE 9.1.

Since the PACS spectrometer continuum data only cover the peak of the spectral energy distribution (SED) for S 140, we have supplemented our analysis with shorter and longer wavelength images. We used the SOFIA/FORCAST images at 11.1, 31.5, and 37 μm discussed by Harvey et al. (2012) and the 24.5 μm Subaru/COMICS data of de Wit et al. (2009). Since all of these observations were obtained with resolution better than that of the PACS/Spec images, we used them both at their native resolution and reconvolved to the PACS/Spec 187 μm resolution. These raw and convolved images were also resampled to a 1′′ grid like the PACS/Spec images. Longward of the PACS wavelengths, the highest resolution image available is in the JCMT SCUBA archive² at 450 μm. We have used this image without further processing in our analysis below. We also examined the SCUBA 850 μm image, which is qualitatively similar but has lower spatial resolution.

2.2. IRAM data

The 30 m observations were conducted during four hours on September 24, 2011. The Eight MIxer Receiver (EMIR) was used during the observations. CO J = 1 → 0 and isotopologues (~110−115 GHz, 3 mm) and CO J = 2 → 1 (230.54 GHz, 1 mm) map data were gathered in two linear polarizations. Our 4′× 4′ maps were centered close to S 140-IRS 1 position at RA = 22:19:18.30, Dec = 63:18:54.2 (J2000) (Fig. 2). For this work, we selected the maps of CO J = 1 → 0 (115.27 GHz), J = 2 → 1 (230.54 GHz), together with the isotopologues ¹³CO J = 1 → 0 (110.20 GHz) and C¹⁸O J = 1 → 0 (109.78 GHz). At all the observed lines, the FTS and WILMA were connected in parallel as backends. In the 3 mm set up, we collected data from three bands of ~ 4 GHz each including the isotopic lines of ¹³CO J = 1 → 0 and C¹⁸O J = 1 → 0 in the CO setup. The beam sizes are ~ 21′′ at 3 mm and ~ 11′′ at 1 mm.

The sky opacity measured by the taumeter at 225 GHz was stable at τ₂₂₅ = − 0.26. Pointing and focus were checked about every hour on quasars and on Uranus and were stable within 2″ and 0.2 mm. The spectral line calibration was checked by pointed observations on IRS 1 and IF. The on–the–fly mapping observations were done in perpendicular directions to avoid scanning artifacts. The typical system temperatures were 96−215 K at 3 mm and 175−190 K at 1 mm. The resulting root mean square (rms) falls between 0.152−0.705 K at 3 mm and ~0.525 K at 1 mm for a velocity resolution of ~0.5 km s^-1 and ~0.25 km s^-1 respectively.

The intensities were converted from antenna temperature units to a scale of main-beam temperature (T_mb), dividing by the main-beam efficiencies η_mb (B_eff/F_eff) of 0.82 at 3 mm and 0.64 at 1 mm. Velocities are given with respect to the LSR. Baselines of orders 0 to 4 have been subtracted. Data reduction and analysis were performed using the GILDAS software³.

2.3. HIFI data

To test the reliability of our method for the gas analysis, we used single-pointing observations (spectra, Fig. 6) and cuts from Herschel-HIFI (Fig. 2, top), in addition to the 4′ × 4′ maps from IRAM-30 m. Single-pointing, frequency-switching spectra in bands 1a, 1b, 2b, 3a, 4a, 4b, 6b, and 7b and OTF maps were analyzed. The single-pointing data were taken toward S 140-IRS 1 at RA = 22:19:18.21, Dec = 63:18:46.9 (J2000) and toward a position in the IF at RA = 22:19:11.53, Dec = 63:17:46.9 (J2000).

For this analysis, we used only the observations of CO J = 9 → 8 (1036.9 GHz), together with the isotopologues, ¹³CO J = 10 → 9 (1101.35 GHz), C¹⁸O J = 9 → 8 (987.56 GHz), toward IRS 1 and IF (obsIds: 1342195050, 1342196426, 1342219209, 1342195049, 1342201741). In addition, we used the cuts of CO 9 → 8 and ¹³CO J = 10 → 9 (obsIds: 1342201741, 1342201806). The observed line parameters in units of antenna temperatures were corrected in units of main beam temperatures using the main beam efficiency of 0.74 (987.56 GHz, 1036.9 GHz, 1101.35 GHz) (Roelfsema et al. 2012), while the beam size at these frequencies is ~ 21′′.

3. Observational results

3.1. Dust continuum observations

For this analysis we selected our 45′′× 45′′ maps in a number of continuum wavelengths with no obvious line emission or with line emission that was easily masked. The particular wavelengths chosen were centered on the wavelengths for which convolution kernels have been developed by Aniano et al. (2011), which enabled us to convolve all maps to a common resolution. These wavelength PACS/Spec data are 73, 75, 84, 94, 110, 125, 136, 145, 150, 168, and 187 μm. We extracted 3% bandwidth continuum fluxes centered on each of these wavelengths while masking out obvious line emission. For comparison with other wavelength data, as well as for comparison with modeling results, we sampled all these images again to a 1′′ spatial grid (Fig. 1). For the final comparison the model images were convolved to a telescope point spread function, and typically a 9′′ square pixel was used to simulate the PACS/Spec observations.

Figure 1 shows the image at 73 μm. The 37 μm (SOFIA) emission is overlaid using black contours showing the positions of IRS 1, 2, and 3. As shown in the image, the emission from IRS 1 dominates at this wavelength.

3.2. The spatial and velocity distribution of CO emission

Figure 2 shows the maps of the peak intensities as observed in their original angular resolution at their peak velocities of CO 1−0, CO 2−1, ¹³CO 1−0, and C¹⁸O 1−0. The lines of the main isotope show two clearly disjointed peaks, one at the densest cluster between the YSOs IRS 1 and IRS 3, and one marking the external interface of the cloud illuminated by HD 211880. Because ¹³CO and C¹⁸O trace the column density structure, the interface peak is shifted into the cloud for these lines, merging with the peak around the embedded cluster.

Finally, while CO 1−0 and CO 2−1 peak toward IRS 1−3 (ΔRA: 10.70′′, ΔDec: –4.0′′), the ¹³CO 1−0 and C¹⁸O 1−0 peak closer to IRS 2 (ΔRA: –11.20′′, ΔDec: 18.30′′). C¹⁸O 1−0 shows a second peak southeast of IRS 1. This reflects a column density effect that is also revealed later in the column density map (see Fig. 11).

Figure 3 shows channel maps of CO 1−0 to illustrate the velocity structure. We find a clear velocity offset between the gas around the cluster and the interface and a much narrower velocity distribution of the interface gas. Figures 4−5 show the line profiles of CO and isotopologues as observed toward IRS 1 and IF. The peak velocity of the emission changes with position on the map. The IF peaks close to the velocity of the source that is ~ −7.1 km s^-1, while IRS 1 peaks at ~ −3.6 km s^-1. The difference may be caused by outflows driven by the infrared sources. For the above reasons we chose to model the peak intensities of the lines in their peak velocities.

3.3. Line profiles

The line profiles as shown in Figs. 4 and 5, are stronger toward IRS 1 and weaker toward the IF (Table 1). At the individual positions, the lines from CO are broader than the lines from the other isotopologues of CO (Figs. 4 and 6).

Table 1 presents the line parameters as measured toward the IRS 1 and IF positions, applying a single-component Gaussian fit, including both HIFI (Fig. 6) and IRAM (Figs. 4 and 5) observations. The uncertainties quoted for the peak intensities correspond to the observational rms, and the uncertainties for the full width at half maximum (FWHM) and υ_LSR are from the Gaussian fits to the peak position. The width of the lines varies throughout the cloud, showing broader profiles toward the center of the map, where the three infrared sources are located and narrow toward the ionization front (Fig. 7). Since we are interested in the quiescent gas and thus the narrow component of the lines, we used observed and computed peak intensities in our calculations, because the narrow component dominates the total peak intensity of the lines (≥50%). With this method we limit the effects from outflow activities on our calculations, but we cannot totally exclude them. The outflow contamination of the line is stronger toward the central sources where the integrated contribution of the broad and narrow components is comparable and weaker in positions away from the sources where the narrow component provides 70−85% of the total peak intensity. We restricted the peak intensity analysis to a velocity window of the width of the FWHM of the lines around the ¹³CO 1−0 line. In this way we always model the same quiescent component.

4. Physical conditions

4.1. RADEX fitting

4.1.1. Method

We used the non-local thermodynamical equilibrium (LTE) radiative transfer program RADEX (van der Tak et al. 2007) to compare the observed line intensity ratios with a grid of models for deriving kinetic temperatures, gas densities, and column densities. As model input we used the molecular data from the LAMDA database (Schöier et al. 2005) and CO collisional rate coefficients from Yang et al. (2010). RADEX predicts line intensities of a molecule for a given set of parameters: kinetic temperature, column density, H₂ density, background temperature, and line width.

In the density range relevant to S 140, the CO 1−0/2−1 ratio is a good tracer for kinetic temperatures. This is demonstrated in Fig. 8, which shows the CO 1−0/2−1 line ratios for a CO column density of 10¹⁸ cm^-2. For typical gas temperatures, we see that these ratios are insensitive to H₂ density for values >10^4.5 cm^-3, well above the critical density of both transitions. At this point we should also point out that at low kinetic temperatures, this ratio is also not a good indicator of kinetic temperature at the accuracy needed to distinguish between 30 K and 60 K. In this temperature range of CO 1−0/2−1 ~ 1.0 ± 0.1, one would have to know the main beam temperature reliably to within 10% to trust the solution to the kinetic temperature to within ±15 K. Given the uncertainties in antenna and receiver calibrations, beam sizes, sidelobes, etc., this level of accuracy is difficult to achieve. One main source of uncertainty is the coupling efficiency of the beam to the source structure. For sources larger than the main beam, the 30 m errorbeams (and sidelobes) contribute to the coupling. At 230 GHz, the main beam efficiency of the 30 m is high, 60%, and the three 30 m error beams contribute 26% to the total power received by the beam pattern. In addition, assuming a reasonable source size of about 1′, only the main beam and first error beam contribute. Because the first error beam contains only 4% of the total beam power, the correction factor for the antenna temperatures, F_eff/B_eff decreases only a little, from 92/59 = 1.56 to 92/63 = 1.46. However, it is not the purpose of our paper to deconvolve the channel maps for the error beams or to provide a full accounting of the error budget contributing to T_mb.

4.1.2. Assumptions

For the CO 1−0, CO 2−1, ¹³CO, and C¹⁸O peak intensities (~ 3 × rms), we performed a χ² minimization to fit the gas kinetic temperatures and column densities assuming that all lines arise from the same gas. For more accurate modeling, the CO 2−1 map was convolved to a lower angular resolution in order to be consistent with the other lines (~21″). The χ² function was computed as the quadratic sum of the differences between the observed and the synthetic line intensities for a range of kinetic temperatures (10 K <T_kin< 200 K) and column densities (10¹³ cm^-2<N_CO< 10¹⁹ cm^-2), which are values that are consistent with the expected ones for such a region (e.g., van der Tak et al. 2000; Spaans & van Dishoeck 1997; Poelman & Spaans 2006; Hüttemeister et al. 1993).

An additional free parameter in RADEX is the temperature of the background radiation field that may pump the line transitions. We adopted the value of 2.73 K for all our calculations since the cosmic background radiation (CMB) peaks at 1.871 mm, and it is the prominent component at millimeter wavelengths. We used a fixed line width value of 3.5 km s^-1 that approximates the value that we have measured throughout the cloud for the narrow component. Finally, to compute the intensities of the isotopic lines, we assume fixed isotopic ratios of ¹²CO /¹³CO = 60 and C¹⁶O/C¹⁸O = 560 (Wannier 1980; Wilson & Rood 1994).

4.1.3. Gas density

Previous studies have shown that the S 140 region contains gas with n_H2~ 10⁵ cm^-3. Poelman & Spaans (2006) report clumpiness in the region and give a value of n_H2 ~ 10⁴ cm^-3 for the interclump medium and ~4 × 10⁵ for the clump gas.

Figure 8 shows that it is impossible to derive the gas density from the low-J lines of CO (2−1, 1−0) for densities of n_H2 ~ 10⁵ cm^-3 or above. A reliable determination of the gas density was thus only possible for those parts of the map where the HIFI cuts provided additional high-J line data. Here, we performed a three-parameter RADEX analysis, fitting the CO 1−0, CO 2−1, CO 9−8, ¹³CO 1−0, ¹³CO 10−9 and C¹⁸O 1−0 lines convolved to the same resolution (21′′), to determine the n_H2. The result is shown in Fig. 9. This fit proves that, at least along the cuts, we find gas densities everywhere, well above n_H2 = 10⁵ cm^-3, where the CO 1−0/CO 2−1 line ratio can be directly used as a temperature measure. In this way we can be sure to avoid the regimes where the ratio is sensitive to the H₂ density.

To assess the reliability of the results, we performed two separate RADEX fits along the cuts where the 9−8 and 10−9 transitions of CO and the isotopologues were observed by HIFI. In the first fitting approach, we fixed n_H2 to 10⁵ cm^-3 and fit only the low-J CO lines observed by IRAM and all the lines (Fig. 10). The second fit uses all the lines and treats the gas density as a free parameter in the range between 7 × 10³ cm^-3 and 7 × 10⁵ cm^-3 as in Fig. 9. The best fit of the latter occurred for n_H2> 10⁵ cm^-3. Adding the higher-J lines drives the fit to systematically higher kinetic temperatures by about ~5−15 K (Fig. 10). The opacities of the lines are found to be in the following ranges: CO 1−0: 16−45, CO 2−1: 52−154, CO 9−8: 3.6−65, ¹³CO 1−0: 0.3−0.8, ¹³CO 10−9: 0.017−0.38, and C¹⁸O 1−0: 0.025−0.09. As an attempt to test how much the optically thin ¹³CO 10−9 line influences the solution, we reran our calculations, applying a double weight for this line. The resulting kinetic temperatures changed by ${\mathrm{~}}_{-1}^{\mathrm{+}\mathrm{2}}$ K, that is, inside the range of the reported errors.

We ran the same kind of analysis toward the two positions with most lines observed⁴ (see Table 1). Toward IRS 1, the full analysis provides a gas kinetic temperature of $\mathrm{~}\mathrm{6}{\mathrm{0}}_{-2}^{\mathrm{+}\mathrm{6}}$ K and a column density of $\mathrm{~}{\mathrm{3}}_{-2}^{\mathrm{+}\mathrm{1}}\mathrm{\times }{\mathrm{10}}^{\mathrm{19}}$ cm^-2, while the same procedure using only the IRAM data results in a kinetic temperature of $\mathrm{~}\mathrm{4}{\mathrm{9}}_{-3}^{\mathrm{+}\mathrm{9}}$ K and a column density of $\mathrm{~}\mathrm{1.}{\mathrm{3}}_{-2.2}^{\mathrm{+}\mathrm{2}}\mathrm{\times }{\mathrm{10}}^{\mathrm{19}}$ cm^-2. Toward the IF position, the full dataset indicates a kinetic temperature of $\mathrm{~}\mathrm{4}{\mathrm{0}}_{-2}^{\mathrm{+}\mathrm{5}}$ K and column density of $\mathrm{~}\mathrm{1.}{\mathrm{0}}_{-1.5}^{\mathrm{+}\mathrm{2.5}}\mathrm{\times }{\mathrm{10}}^{\mathrm{18}}$ cm^-2, while the IRAM data only result in a kinetic temperature of $\mathrm{~}\mathrm{3}{\mathrm{5}}_{-8}^{\mathrm{+}\mathrm{9}}$ K and column density of $\mathrm{~}\mathrm{4.}{\mathrm{4}}_{-2.5}^{\mathrm{+}\mathrm{3.6}}\mathrm{\times }{\mathrm{10}}^{\mathrm{17}}$ cm^-2. The gas temperatures obtained when using only the low-J lines from IRAM underestimate T_gas, providing a lower limit of gas temperatures (Yıldız et al. 2013). The major and minor axes of the resulting χ² contours were used as the error bars of the two parameters.

For a more accurate determination of the H₂ density, tracers such as CS 2−1/3−2 and HCO⁺ 1−0/3−2, are more reliable (van der Tak et al. 2007), but were not observed. Snell et al. (1984) derived a density of 5 × 10⁵ cm^-3 toward IRS 1−3 using four CS transitions that agree with the value we derived. In addition Goldsmith & Langer (1999) performed a LTE population diagram analysis toward IRS 1 using the same CS dataset and derived a kinetic temperature between 30−50 K, which is lower than our value probably due to the effect of their larger beam size.

4.1.4. Kinetic temperature and column density distribution

The resulting maps of kinetic temperatures and CO column densities are shown in Fig. 11. We find a kinetic temperature of ~ 55 K toward the center and ~ 45 K toward the ionization front, while the rest of the cloud is characterized by lower temperatures (25−40 K; Fig. 11). The column density toward the center was found to be the highest, with a value of ~ 6.5 × 10¹⁸ cm^-2, while toward the IF was found to be ~ 1.1 × 10¹⁸ cm^-2. The lowest column density that was determined throughout the cloud is ~ 1.4 × 10¹⁶ cm^-2. The gas temperature map obviously reflects the different heating contributions by the embedded cluster of the three high–mass YSOs (IRS 1 to 3) in the center and from the external B0V star HD 211880 toward the IF.

4.2. Dust analysis

Prior to any radiative transfer modeling of the dust in this region (Sect. 5), we here estimate the rough temperatures and optical depths relevant to the dust distribution. For this we used the PACS/Spec continuum data from 73−187 μm, which covers the peak of the SED over the entire region mapped with PACS well (Fig. 12). To analyze our dust continuum observations, we used two subsets of the PACS and SOFIA images. To compute a luminosity map, we used the 11−187 μm images all as reconvolved to the 187 μm resolution (13′′). At each 1′′ pixel of the computed image we integrated the flux density from 11.1−187 μm and included a linear extrapolation to zero flux level beyond 187 μm, which typically results in only a few percentage points added to the total. We believe the absolute uncertainties in this map are approximately 15% based on the absolute calibration uncertainties of all the input data. The relative uncertainties are probably much less, ±5%, since there is very little luminosity shortward of 10 μm and longward of 200 μm.

We computed color-temperature and optical-depth values at each spatial point on the 1′′ grid by fitting a blackbody function modified by a λ^-1 dust emissivity variation. The peak dust optical depths at the shortest wavelength are likely to be several × 0.1, but still less than unity, so we did not include any effects from optical thickness in these estimates. A comparison of the luminosity, color temperature, and optical depth maps and several images is shown in Fig. 12. These images show a number of important qualitative properties of the dust emission. First, there is a clear and relatively smooth change in source morphology from that at 73 μm, which is similar to the 10−37 μm images, to the morphology at 187 μm, which is beginning to show many of the features of the SCUBA 450 μm map. Second, the dust optical depth we derived from the Herschel data is very similar to the 450 μm emission map. Finally, the luminosity image shows clearly that IRS 1 is the main luminosity source in the region, followed in importance by IRS 2 and then IRS 3. There is no obvious luminosity peak within the area of strongest 450 μm emission, so this is probably a peak in the column densities of dust and gas.

One uncertainty of the dust temperature may stem from the assumed spectral properties of the dust grains, in particular the assumed spectral index β. Recent observations of dust in the diffuse ISM with the Planck mission indicate a mean dust temperature of ~20 K and dust emissivity index of ~1.6 (Jones 2014). This is clearly higher than the value of β = 1 assumed here, which is more appropriate for dense clouds (Ossenkopf & Henning 1994). Since our region is denser than the majority of those observed with Planck, we can consider the value of β = 1.6 as an extreme upper limit for our case. When we fit a blackbody function modified by a λ^-1.6 to derive the dust temperatures, we obtain temperatures that are lower by 5−20 K than the values from the adopted λ^-1. We conclude that in our analysis we derive an upper limit of dust temperatures.

4.3. Comparison of dust and gas temperatures.

For a direct comparison with the gas kinetic temperature map, the dust temperature map was convolved to the same angular resolution (~21″). The comparison between gas and dust temperatures from the low-J line analysis (Figs. 10 and 13), shows a very similar spatial structure around the infrared sources. The higher-J lines from HIFI, though, reveal that the low-J analysis underestimates the gas temperatures (see Sect. 4.1.3). Assuming that this trend applies to the entire region and not only along the cuts (Fig. 10), we conclude that the gas temperature is higher than the dust in the entire region.

The temperature difference of ~5−15 K revealed from the complete analysis lies outside the uncertainty range of the 2−7 K and indicates a more efficient gas heating even at densities ≥10^4.5 cm^-3 where the two components are relatively well coupled.

To illustrate the expected overall temperature dependence, we show in Fig. 14 the outcome of a simplified PDR model for S 140 computed with the KOSMA-τ PDR code (Röllig et al. 2006, 2013). It shows the temperature profile for a spherical PDR clump with a mass of 100 M_⊙, mimicking a plane–parallel PDR, and a surface density of 10⁵ cm^-3 illuminated by an UV radiation field of 100 Draine fields (Draine 1978). For reference, the figure shows the visual extinction A_V as a function of the depth into the cloud. The model shows that through almost all of the whole cloud, the gas temperature is higher than the dust temperature by about 15 K for A_V = 0.1 but by only 2 K deep in the cloud. The small kink in the gas temperature around A_V = 1 stems from H₂ formation heating on PAHs and the temperature structure of the PAHs. Overall, the model confirms that for a homogeneous medium, we expect the same temperatures for gas and dust within our observational measurement errors for visual extinctions into the cloud of more than about 0.5, i.e. for depths of one arcsecond and more at the distance of S 140. The extended hot gas therefore cannot be explained by a homogeneous extinction of the UV radiation but requires some clumpy structure, allowing for deeper UV penetration.

5. Effect of density gradients

Assuming uniformly distributed densities and temperatures for each spatial position was a good approach for the extended gas and dust in S 140, but detailed modeling is required for the more complicated dust and gas structures that are expected to be close to the infrared sources. The radiation we receive is a result of various processes, including the interaction with the surrounding dust. For a more precise physical approach, the radiative transfer calculations should include the original radiation, characteristics of the dust grains and the dust density distribution. Density gradients have been revealed for S 140 in previous studies that include Harvey et al. (1978), Gürtler et al. (1991), van der Tak et al. (2000), and Maud & Hoare (2013). We first performed an advanced dust modeling taking dust density gradients into account. Then we applied the best-fit derived dust mode to predict the CO intensities using an advanced radiative transfer code (i.e., RATRAN; Hogerheijde & van der Tak 2000), as an attempt to test the accuracy of our dust model comparing with the gas observations and predictions.

5.1. Dust modeling: an approach

It is clear from all the infrared and sub-mm images of the S 140 cluster that there are several luminosity sources in the central arcminute and that the dust distribution is unlikely to be spherically symmetric. This implies that any simple model of the dust heating will be limited in its applicability. We have, however, identified some goals to address with simple radiative transfer models. First, since all the observational data suggest that IRS 1 is the dominant luminosity source, we believe that we can determine some rough properties of the dust distribution close to it by ignoring the heating from IRS 2 and 3 and any of the other much lower luminosity objects that are part of the embedded stellar cluster around IRS 1. To this end, we concentrate on two observables for the model comparison for IRS 1: (1) the SED of IRS 1; and (2) the spatial flux distribution from the center to the south of IRS 1 where there are no strong obvious nearby heating sources seen in the infrared that are likely to contribute significantly to the dust heating. A second goal is to try to understand whether the peaks in the 450 μm image can be due to some unusual dust distribution to the west and southwest of IRS 1 with no additional internal heat source or if some internal heating is required to explain them. For example, the recent observations and modeling by Maud & Hoare (2013) suggest an internal source within the strongest 450 μm peak to the southwest of IRS 1, i.e. SMM 1 (ΔRA: –14.22′′, ΔDec: –7.90′′, relative to IRS 1).

The approach we took to addressing these goals involved several steps that we discuss in detail in the following sections. First, we tried to find what properties of models were required to provide a fit to the observed SED from the central pixel centered on IRS 1 and to the source profile of IRS 1 to the south, assuming a single, spherically symmetric dust distribution. This model was then used as a starting point for subsequent stages. The second step was to create a model made from hemispheres of two different spherically symmetric dust distributions. Though this model is simplified, it seems likely that far from the boundary region, such a combination can give insight into the degree to which such a dust distribution might reproduce the observed source structure to the west and southwest of IRS 1 by assuming a higher column density in this region than to the south and east of IRS 1. For this comparison we used the source profile extending southwest along a line from IRS 1 to and beyond the position of SMM 1, as given by Maud & Hoare (2013). As an alternative to the two-hemisphere model, we also attempted to fit the profiles along this line with two separate sources by using the superposition of two spherically symmetric models with different central source luminosities and very different dust distributions separated on the sky by ~16′′ to attempt to reproduce the observed 1D flux profile between IRS 1 and the position of SMM 1 southwest of it.

In testing these last two model variations, the two-hemisphere construction and the two-source construction, it became clear that the best results would probably result from a combination of both features. Such a model would also be most realistic in light of the much higher dust column density to the west and southwest as indicated both by the SCUBA maps and our Herschel maps, together with the presence of the compact sub-mm/mm source found by Maud & Hoare (2013). After constructing a few such models, we realized that with so many free parameters, it would be difficult to find the range of well-fitting parameters without testing the fits over a large, multi-dimensional grid. This process represented the final stage of modeling for IRS 1, and we now describe the process and results. Figure 15 shows a schematic diagram of the model for IRS 1 and the sub-mm peak to the southwest.

For all the models we used the DUSTY code (Ivezic et al. 1997) and converted the dimensionless output values to those appropriate for the assumed distance of S 140, 746 pc (Hirota et al. 2008) and the assumed luminosity (discussed below) for comparison with the observations. We tried models with two different dust compositions: the first was a mixture of 90% DL silicates (Draine & Lee 1984) and 10% amorphous carbon, both of whose optical constants are distributed with DUSTY. The second was OH5 (Ossenkopf & Henning 1994) dust, for which a number of authors (e.g., van der Tak et al. 1999) have suggested that they may be more appropriate in dense star-forming regions where significant grain growth may have taken place. We found similar quality fits with comparable dust density distributions for both grain compositions. The 1 μm optical depths were, however, about twice as high with the DL/carbon dust than with the OH5 grain properties owing to the overall difference in the slope of the dust optical depth between 1 and 100 μm. The output of the “observed” profile of each model at each wavelength was first interpolated onto a 3050 × 850 grid with spacing equivalent to 0.05′′. We then convolved this image with the appropriate observing beam to compare with the observed images (not those convolved to the 187 μm resolution). Finally, we “observed” this image with a 9′′ square beam that we moved across the image in 1′′ steps to compare with the observed images whose fluxes we also added into such a moving 9′′ square beam. For the shorter wavelength images, we also applied the same 9′′ square pixel flux summation for comparison between the models and observations.

5.2. Model grid

The model that we used is described by nine parameters, four for the eastern half of IRS 1, four for the western half, and one, luminosity, for SMM 1. Table 2 lists these nine parameters, as well as the range of values tested in a grid, thereby allowing all possible combinations of these values (150 000 models). A preliminary run of ~800 000 models showed that the results were quite insensitive to the value of parameter 4, the radius where the slope of the density gradient changed in the high-column hemisphere, so it was fixed at a value of 300× the inner radius. Similarly, the value of the density gradient in this inner region affected the final results only slightly, so this was fixed at a value 0.4 (ρ ∝ r^-0.4). Therefore, for the final results there are seven free parameters.

To select the most likely models and to estimate the range of values that produces a reasonable fit, we characterized the model fits by the χ² that we computed as follows. We parameterized the observations and accompanying model fits into 30 values of the SED and the source profile at various wavelengths. Since the signal-to-noise ratio of almost all the observations is quite high, the systematic uncertainties of the observations (e.g., pointing) dominate the true uncertainties. We, therefore, assigned somewhat arbitrary values to the assumed errors for these 30 values to drive the fitting to something “reasonably” close to the observations. Table 3 lists these properties of the observed data and uncertainties.

As shown in the table, the parameters used to calculate the χ² of the fit included the SED between 37 and 450 μm, the flux at the position of SMM 1, the relative flux at 37−450 μm at offsets of 9′′ and 18′′ to either side of the position of IRS 1 (i.e., one pixel and two pixels), and finally, the observed luminosity at those four offset positions (not relative to the central pixel). We intentionally did not include the portion of the SED shortward of 37 μm, since the dust close to the central source of IRS 1 is likely to be distributed in a disk-like configuration, which would lead to more flux escaping at shorter wavelengths than is consistent with the assumption of spherical symmetry. Likewise, we did not attempt to fit the total observed luminosity at the central pixel since much of the luminosity is emitted shortward of 37 μm. To estimate the range of model parameters that provided a “good” fit to the data, we assigned probabilities to each model equal to exp( − χ²/ 2).

The best-defined parameter values from the probability plots are the mild density gradient away from IRS 1 of ρ ∝ r^-0.4 and the outer radius of 1500 times the inner radius. The optical depth is relatively well defined by A_v ~ 40, though there is a better-defined joint probability that includes the optical depth and density gradient. The peak probability for the luminosity of SMM 1 is 100 L_⊙, but there is a wide range of values with reasonable probability from ~50 L_⊙ to 300 L_⊙. The parameters of the dense material to the southwest of IRS 1 are not well defined for the most part. Interestingly, the gradient of the material in the outer area (300−1500 × the inner radius) seems to have a very probable range similar to the range found for the rest of the cloud around IRS 1; i.e., r^{− p} with p ~ 0 − 0.5. There is a significant joint probability distribution between the outer gradient to the southwest and the luminosity of SMM 1. This is to be expected, since the effects of a steeper gradient can be largely compensated for by increasing the luminosity of SMM 1. Although the most probable models have A_v = 40, the probability distribution suggests that the most probable range is 40−60, i.e., somewhat higher than the gradient to the south.

Examples of the model-fit profiles from one of the two lowest χ² (2.7) fits, as well as the SED, are shown in Figs. 16−18. The model profiles shown in Fig. 16 display a less than perfect fit at the longest wavelengths in our PACS observations, 140−190 μm. This effect has been driven by two parts of the fitting process: that we did not include the relative flux at the central (IRS 1) position in the χ², and that it proved quite difficult to find any models that dropped as slowly as the observations at the ends of the profiles 18′′ away from the center. The properties of the models that do the best job fitting the profiles 18′′ away from IRS 1 also drive the relative 187 μm profile to be too high at IRS 1 compared to its value at SMM 1. In some simple tests, we have found that if we added an artificial floor to the model fluxes beyond 125 μm, we could fit the longer wavelength profiles significantly better.

5.3. IRS 2 and 3

The most important factor affecting our attempts to model IRS 2 and 3 is that the flux from our well-fitting models for IRS 1 is a significant fraction of the total flux at the positions of IRS 2 and 3. This means that any models we make for IRS 2 and 3 will have an additional large uncertainty due to the uncertainty in the true contribution from IRS 1. For example, the model flux from IRS 1 at the position of IRS 3 is more than 50% of the total at wavelengths beyond 120 μm. We have therefore computed models for roughly a dozen possible configurations for each source to develop some feeling for the range of likely parameters, but the uncertainties in our estimates are quite large.

For IRS 3 we used a profile that begins on the east side of the source and, after passing through the center, goes straight to the south in order to minimize the confusion from IRS 1. An additional complication is that IRS 2 is clearly elongated in the east-west direction in the mid-IR, and there is a clear position shift in the same direction between the mid-IR and our PACS data. This suggests that there are at least two relatively luminous objects heating the dust at the position of IRS 2 with different dust optical depths, higher to the east and lower to the west. The position shift between 12 and 73 μm is only a few arcsec, which is a fraction of our PACS native pixel size. Therefore, for IRS 2 we used a source profile from east to west and a model like the one in the “two-hemisphere” modeling for IRS 1, with a hot, lower-optical-depth dust distribution on the western side and a cold, higher-optical-depth distribution offset 4′′ to the east. Such a model will not accurately reproduce the source elongation, but should roughly fit the overall SED and some part of the profiles. Figures 19−22 show the model SED and profile fits for some representative models that illustrate both the features that can be fit reasonably well as for those that are difficult to fit.

For both of these sources, we did attempt to fit the mid-infrared fluxes in the SED, unlike for IRS 1. The very extended tail of emission on the west side of IRS 2 at 125 and 187 μm coincides well with the increase in optical depth seen in our data and the SCUBA maps in that area. We have not tried to fit this with these simple models. The bump seen on the eastern side of the long wavelength profiles for IRS 3 is coincident with the outermost pixel of our mapping, so we did not try to fit it.

The most well-defined results seem to be that 1) the dust distribution around both sources is likely to be quite flat; and 2) the luminosities of the sources are likely in the range of 1000−2000 L_⊙, as shown in Table 4. The first could be consistent with very little dust around them that is associated with them, but mostly with the extended distribution around IRS 1.

5.4. Gas modeling

To test our best-fit DUSTY model, we ran the Monte Carlo radiative transfer code RATRAN (Hogerheijde & van der Tak 2000), which treats the physical structure of the sources, including temperature and density gradients. RATRAN estimates the local radiation field at all line frequencies taking the radiation field from every other position in the cloud into account. We ran RATRAN applying the best-fit DUSTY model and assuming the same temperature for dust and gas toward IRS 1 in order to model the CO 1−0, CO 2−1, CO 9−8, and isotopologue lines and to compare them to the outflow–subtracted component from the observations (Fig. 23). The outflow emission was subtracted since the DUSTY models focus on the bulk of the quiescent dense gas where the outflow cavities are negligible. For the outflow subtraction, we applied a two-component Gaussian fit to the observed lines toward IRS 1 followed by subtraction of the broad component.

For our calculations we defined a grid of 18 spherical shells for the eastern (lower density) hemisphere (constant ρ ∝ r^-0.4) and 22 for the western (high density) hemisphere by applying the density gradients as obtained from the DUSTY model. The inner and outer radii were taken from the DUSTY model and were set to 2.08 × 10¹⁴ cm and 3.12 × 10¹⁷ cm, respectively, while the inner H₂ density was set to 9 × 10⁵ cm^-3. This value was derived from the best-fit dust model assuming that the gas is entirely molecular and using a mean gas mass per hydrogen of 1.4 amu and a gas-to-dust ratio of 100.

Dust continuum radiation is considered using the same OH5 opacities as in the DUSTY model. We assume a static envelope without infall or expansion, a fixed CO abundance of 10^-4, and a turbulent line width of 3.5 km s^-1.

Figure 23 shows the lines of CO 1−0, CO 2−1, and CO 9−8 and isotopologues as observed toward IRS 1 overplotted with the convolved (~21′′) synthetic emission as calculated with RATRAN. We show the resulting line profiles for both the low-density east hemisphere and the high-density west hemisphere from the DUSTY models. With the exception of the low-J transitions of the isotopologues, we observe no significant differences between the two cases. The fit between observed and modeled lines shows that the low-J lines are reproduced but the higher-J lines are significantly underestimated, especially for the isotopologues. The model cannot reproduce the highly excited lines tracing high column densities of warm and very dense gas. An adapted kinematic structure and/or clumpiness can potentially remove the modeled self-absorption dip of the optically thick ¹²CO lines, but treating these effects is beyond the scope of this work. Models with T_gas>T_dust, as inferred from the simple analysis in Sect. 4.3, would probably be able to reproduce the high-J lines without significantly affecting the low-J lines. Local density enhancements can be an alternative solution, if the total mass of the cloud is conserved, such as in a clumpy or disk/outflow geometry.

6. Discussion and conclusions

Based on PACS, HIFI, and IRAM data, we find a wide range of dust and gas temperatures in the S 140 region. The warmest gas (~ 55 K) is around the three infrared sources, while the surrounding environment is characterized by colder temperatures (~25−35 K). The gas temperature increases toward the southwestern edge of the cloud (IF), reaching values of ~40−45 K. This rise in temperature is the result of the UV heating from the external B0V star HD 211880 that lies ~7′ southwest of the edge of the cloud. Unfortunately, we cannot compare the gas and the dust temperatures toward IF because of the smaller area that our continuum dust observations cover.

We find that the dust density gradient around IRS 1 is likely to be shallower than the best fits found in some earlier models in order to explain the spatial structure we find in the far-IR. The strong emission to the west and southwest of IRS 1 at wavelengths longward of 100 μm can possibly be explained by heating from IRS 1, along with a strongly increasing density to the west and southwest. We do, however, find a significantly better fit by interpreting the emission as arising from two or more separate sources with very different amounts of dust around the sub-mm source compared to IRS 1. The luminosity of the internal source at the sub-mm peak is likely to be a few ×10²L_⊙ in this model. This model with an internal luminosity source at the sub-mm peak is consistent with the study of Maud & Hoare (2013), who find that the two strongest 1.3 and 2.7 mm emission peaks are the position of IRS 1 and SMM 1. Because our models did not include outflow cavities or clumpiness, it is possible that this conclusion should be revisited after more extensive modeling.

The gas temperature analysis, including both high and low-J lines, was possible only along the HIFI cut and revealed a systematic excess of gas temperatures against dust temperatures. This result indicates a more efficient gas heating even at densities ≥10^4.5 cm^-3 where the gas and dust are usually expected to be in thermal equilibrium. New SOFIA/GREAT mapping observations of CO 13−12 and 16−15 confirm this excess of the gas temperature, and the results show an even greater gap in energy between the gas and dust components (Ossenkopf et al. 2015). The high-J lines show that a limitation to the low-J analysis (IRAM) causes an underestimate of the gas temperatures by ~7−12 K and thus provides only a lower limit of the temperatures in the cloud. Thus it is most likely that the gas temperature is higher than the dust also in the whole field. Unfortunately, the higher-J lines do not cover the entire region so that we cannot prove this. The detailed gas modeling also indicates that DUSTY results are good estimates of the temperature and density structure around IRS 1 for the colder gas, but there is an excess of warm dense gas that cannot be reproduced by these models (i.e., CO 9−8). The attempt to fit those higher transition lines by modeling hotter gas than dust and/or applying different geometries (i.e., clumps, disk geometry) that would increase the density locally but not the total amount of mass is in our future plans.

The observed gas temperature excess cannot be explained by a higher cosmic ray ionization. Van der Tak & van Dishoeck (2000) reported a cosmic ray ionization rate of <10^-16 s^-1, which is too low to provide a major gas heating in S 140. Some gas heating may be due to outflow activities. In the RADEX and RATRAN analysis, we focused on quiescent gas and limited the contribution from protostellar outflows. However, it is possible to have high-velocity motions perpendicular to the line of sight. In addition, the larger difference between gas and dust appears in points close to the observed outflow by Preibisch & Smith (2002). A collection of shocks, longitudinal and transversal to the walls of the possible inner cavity around the IR sources, will contribute to the emission at the ambient velocity. Oblique shock components are thus a plausible scenario for a somewhat more efficient heating of the gas than of the dust.

The most probable scenario appears to be the deep UV radiation in a clumpy medium. Dust is heated by UV and IR, while gas is heated by processes driven only by UV with the main one being the photo-electric heating (Hollenbach & Tielens 1997; Röllig et al. 2013). In the vicinity of any stellar or protostellar source, the gas is much hotter than the dust when dust and gas are not coupled. Dust-gas collisions try to equilibrate them if the density is high enough. The dust quickly attenuates the UV so that the gas becomes colder when going away from IRS 1. The penetration of the IR is deeper (lower extinction in IR compared to UV), so that the decay in the temperature of the dust is shallower than for the gas. However, since we observe that the gas remains hotter than the dust over a relatively large distance from IRS 1, the UV seems to penetrate much more deeply into the cloud than expected from a homogeneous medium. The natural explanation is a clumpy medium where there are always enough rays between the clumps that still allow UV propagation to keep the gas warm.

Previous studies already report that S 140 is clumpy (e.g., Kramer et al. 1998). Battersby et al. (2014) performed a similar study to determine the relationship between gas and dust in a massive star-forming region (G 32.02+0.05) by comparing the physical properties derived from each, using NH₃ transitions. In that study similar temperature differences are reported with the dust temperatures being lower than gas temperatures (by a few K) in the quiescent region, indicating that gas and dust might be not well-coupled in such environments. Differences between gas and dust temperatures in several environments with high densities have been reported in more studies (Papadopoulos et al. 2011; Hollenbach & Tielens 1999).

Our observations trace warm and cold dust but also warm and less warm gas, because we have multiple transitions of CO. The low-J lines tend to show the surface of clumps since they are optically thick but the optically thin isotopologues ¹³CO 1−0 and C¹⁸O 1−0 carry information from deeper positions. Furthermore, the higher transitions CO 9−8 and ¹³CO 10−9 trace warmer gas, but they also arise from deeper positions. That we have optically thin lines and the dust emission is optically thin confirms that we trace the whole line of sight both in gas and dust and thus the observational limitations can probably be excluded.

---

Acknowledgments

The authors would like to thank the referee for a careful reading of the manuscript and for providing helpful comments and suggestions that improved the paper. P.H. was supported through the NASA Herschel Science Center data analysis funding program, by a NASA contract issued by the Jet Propulsion Laboratory, California Institute of Technology to the University of Texas. V.O. was supported through the Collaborative Research Center 956, subproject C1, funded by the Deutsche Forschungsgemeinschaft (DFG). A.F. thanks the Spanish MINECO for funding support from grants CSD2009−00038 and AYA2012−32032. 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 Astrobiologa (CSIC-INTA). Sweden: Chalmers University of Technology – MC2, RSS and GARD; Onsala Space Observatory; Swedish National Space Board, Stockholm University – Stockholm Observatory; Switzerland: ETH Zurich, FHNW; USA: Caltech, JPL, NHSC.

References

All Tables

All Figures