Free Access
Issue
A&A
Volume 563, March 2014
Article Number A97
Number of page(s) 31
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201322541
Published online 17 March 2014

© ESO, 2014

1. Introduction

The chemical evolution of star-forming regions is an important topic from various perspectives. From an astrochemist’s point of view, one wishes to understand the chemical properties of the gas to be able to characterize the chemical composition and evolution of the interstellar medium (ISM) in general. This starts with the chemistry associated with diffuse clouds and the formation of relatively simple molecules like H2 or CO, and continues to the formation of complex molecules (even pre-biotic ones) in hot molecular cores or protostellar accretion disks. Complementing this approach, astrophysicists need to have a good understanding of the chemical constituents of the ISM to use different molecules as tools to probe physical conditions such as densities, temperatures, or kinematical properties. These two approaches are interdependent.

Chemical inventories of (high-mass) star-forming regions have so far largely concentrated either on complete line-surveys toward selected regions (e.g., Orion-KL or NGC 6334I &I(N), Sutton et al. 1985; Schilke et al. 1997, 2001; Walsh et al. 2010) or have targeted small spectroscopic setups toward specific subsamples (e.g., hot molecular cores, HMCs, or infrared dark clouds, IRDCs, Hatchell et al. 1998; Bisschop et al. 2007; Vasyunina et al. 2011). In addition to these line surveys, a few interferometric high-spatial resolution studies of limited samples of high-mass star-forming regions exist (e.g., Blake et al. 1996; Beuther et al. 2009).

Many studies suffer either from too limited sample sizes or spectral coverages to be able to characterize the chemistry of high-mass star-forming regions at various evolutionary stages in a statistical sense. A number of more recent studies aimed at a deeper understanding of a broader chemistry in the various evolutionary stages of high-mass star formation. Vasyunina et al. (2011) compared abundances in the IRDCs with the data for more evolved sources taken from the literature, and the main problem was systematic, where different spectral setups, beam sizes, and integration times were used in different samples. Fontani et al. (2005) studied a larger sample of protostellar candidates from the IRAS Point Source Catalog on the basis of their IR color and observed them in CS and C17O and the 1.2 mm continuum. Based on these observations and results from the literature, they were able to classify the evolutionary stages as very early and probably previous to the formation of an UCHii region. Zinchenko et al. (2009) observed five high-mass star-forming regions in a handful of simple molecular species to estimate their physical and molecular parameters. They found systematic differences in the distributions and the abundances of various molecules and discussed the importance of HCO+, HNC and, especially, N2H+ as potentially valuable indicators of massive protostars. Pirogov et al. (2007) mapped twelve high-mass star-forming regions in CS, N2H+ and 1.2 mm dust continuum and derived the physical parameters, the density and chemical structure of the associated dense cores within. Reiter et al. (2011) also studied 27 high-mass clumps in a larger set of molecules and amongst others discussed the dependence of various physical clump properties on chemical properties. They found that molecular column densities are only weakly correlated with any of the derived physical properties of the clumps. Another still ongoing study is the MALT90 survey, which aims to chemically characterize dense molecular clumps (Foster et al. 2011; Sanhueza et al. 2012; Hoq et al. 2013; Jackson et al. 2013).

To continue the chemical characterization of high-mass star-forming regions, it is important to study all evolutionary stages in a systematic and consistent way for all sources. Furthermore, it is desirable to include more molecular species as well as be able to understand the complete chemistry along the evolutionary sequence with models. In this work we divide the early evolution of the high-mass star-forming regions into four different stages that are observationally motivated and guided by the evolutionary sequence shown in Beuther et al. (2007a) and also in Zinnecker & Yorke (2007) who divided the different stages based on their physical conditions. In the first, earliest considered phase, quiescent infrared dark clouds (IRDCs) are formed. They consist of cold and dense gas and dust and emit mainly at (sub-)millimeter wavelengths. Their physical conditions are close to isothermal. In our sample this group consists of starless IRDCs as well as IRDCs already starting to harbor point sources at μm-wavelengths.

Previous theoretical works, for instance by Krumholz & Thompson (2007); Narayanan et al. (2008); Heitsch et al. (2008), suggested that there might be a long-lived pre-IRDC massive clump stage. This stage was also seen in the observations by Barnes et al. (2011) and supported by extragalactic works, for example, by Koda et al. (2011) who suggested a long-lived pre-MSF stage in M51. The existence of a long-lived pre-IRDC phase as a precursor is not doubted in our approach. But from a modeling point of view, in our opinion ~20 000 years are enough to convert almost all initially atomic gas in regions with densities 105 cm-3 into molecular-rich gas and thus reach the molecular IRDC stage. Thus, based on chemistry, we define the year zero in our evolutionary sequence when the objects reach densities of 104 cm-3 and become detectable as cold dense molecular clouds. The evolution prior to this stage is not considered in this work and took already an unknown amount of time.

In the second phase, the so-called high-mass protostellar objects (HMPOs) form, hosting an actively accreting protostar(s) with >8M, which already shows an internal emission source(s) at mid-infrared wavelengths. The central temperature rises and the temperature profile starts to deviate from the isothermal case. In this work we split the physically homogeneous HMPO phase into two subsequent phases with the HMPO being followed by the much warmer hot molecular core phase (HMC). This phase is distinguished from the earlier HMPO phase from a chemical point of view. In the HMC stage the central source(s) heats the surrounding environment, evaporating molecular-rich ices and giving rise to molecular complexity in the gas. Finally, the UV-radiation from the embedded protostar(s) ionizes the surrounding gas and an ultra-compact Hii (UCHii) region is formed (fourth stage). In these objects many of the previously formed complex molecules are no longer detected because they are most likely destroyed by the ionizing radiation. Physically, many UCHii regions have probably stopped accretion at this point, but this class is not entirely homogeneous either and accretion may still continue in some UCHii regions.

Because the evolution in high-mass star formation takes place on rather short timescales and in clustered environments, the transitions from one into the next stage are smooth and not always clearly distinct. There might be overlaps among the HMPO, HMC and even UCHii region stage. Some regions are already surrounded by an ionized medium but still show the rich hot core chemistry. Due to the large beam, regions with different evolutionary stages in their proximity can also be observed as one object showing the characteristics of the different stages. In this work we wish to disentangle this from an observational point of view to be able to characterize the evolution in every considered stage statistically and also to model the evolution on a consistent timeline. On the one hand, we follow the evolution based on physical quantities, especially the temperature, which is rising from IRDCs to HMPOs to UCHii regions and on the other hand, based on the chemistry which splits the HMPOs into early, chemical poorer HMPOs, and HMCs with a molecular richer chemistry.

To set the observational results into context, an extensive multidimensional modeling and fitting of the observationally derived column densities is performed. The modeling of each individual evolutionary stage is based on the iterative fitting of the data with a set of 1D power-law physical models coupled to a pseudo-time-dependent gas-grain chemical model, calculated over Myr. The cloud density and temperature structures of the environment, as well as its chemical age, are the variable parameters. The chemical characterization of the presumed evolutionary sequence in a statistical sense, together with the unique approach to fit the observed chemical evolution in high-mass star formation directly with a model, enables us to set the results into a broader context.

We describe the source sample in Sect. 2 followed by the observations in Sect. 3. The derivation of the chemical abundances and a chemical characterization of the four different subsamples based on the observational data are given in Sect. 4. In Sect. 5 we describe our model and the fitting process and interpret the observed values in combination with the fitting results in more detail and compare the results with other studies. We conclude with a summary in Sect. 6.

2. Source sample

The source sample contains 59 high-mass star-forming regions, consisting of 19 IRDCs and 20 HMPOs as well as 11 HMCs and 9 UCHiis. The complete sample is listed in Table 1 with coordinates and distances. The telescope pointings were centered on the column density peaks derived from the dust continuum observations obtained either with Mambo (at 1.2 mm) or from the ATLASGAL and SCUBA surveys at 870 μm and 850 μm, respectively (Schuller et al. 2009; Di Francesco et al. 2008). From these maps we determined the corresponding H2 column densities (see Sect. 4.3). The sources are mostly located within the galactic plane, with an average heliocentric distance of ~ kpc.

The sources were selected from different source lists. The lists of the IRDCs were first presented in Carey et al. (2000) and Sridharan et al. (2005) and are part of the Herschel guaranteed time key project EPOS (Early Phase of Star Formation, Ragan et al. 2012). This sample consists of 6 IRDCs showing no internal point sources shortward of 70 μm and 13 IRDCs that have internal point sources at 24 μm and 70 μm (see source list in Table 1 for individual internal point source detections). The HMPOs were taken from the well-studied sample by Sridharan et al. (2002), and Beuther et al. (2002a,b). HMC sources are selected from the line-rich sample of Hatchell et al. (1998) including a few additional well-known HMCs, W3IRS5 and W3(H2O) and Orion-KL. For the UCHiis, we selected line-poor high-mass star-forming regions from Hatchell et al. (1998), and additional sources from Wood & Churchwell (1989).

3. Observations with the IRAM 30 m

3.1. Data collection

The 59 sources were partly observed in the winter semester (20.1.2011-17.3.2011) and mainly in the summer semester (28.10.2011-30.10.2011, 23.11.2011) with the IRAM 30 m telescope at Pico Veleta (Spain). We performed a spectral line survey in two setups at 3 mm and 1 mm using the EMIR receiver with the FFTS backends. The system temperatures of the observations taken in the two different semesters are comparable with Tsys ~ 100 K in the 3 mm band and Tsys ~ 250 K in the 1 mm band. The large 16 GHz bandwidth allowed us to efficiently study the chemical properties of this relatively large sample of high-mass star-forming regions over a broad range of simultaneously observed molecular transitions. Due to changes in the instruments and thus the coverage in frequency during distinct observation runs we only analyzed the frequency ranges at which all sources were observed, namely 8694 GHz, 217221 GHz, and 241245 GHz. The frequency ranges were chosen to cover transitions of important molecules containing nitrogen, oxygen, sulfur, carbon, and silicon. The observations were carried out in wobbler-switching mode with 1.25 min on-source integration time with a spectral resolution of ~0.3 km s-1 at 1 mm and ~0.6 km s-1 at 3 mm. The beam sizes of the IRAM 30 m telescope are 11″ at 1 mm and 29″ at 3 mm. Typical 1-sigma rms values are on the order of ~0.1 K at 1 mm and ~0.03 K at 3 mm, respectively. To reduce all spectroscopic data, the standard GILDAS1 software package CLASS was used.

3.2. Data analysis strategies

For the analysis we focused on a specific sample of 15 different molecular species for which we calculated the column densities (see Table 2 for a list of the molecules). To directly compare the derived values with the modeled values (see Sect. 5.1) we relied on the most common molecules among the different sources at a distinct evolutionary stage (e.g., N2H+, HCN, C2H, HNC, HCO+). These molecules we complemented with more complex molecules that are typical of more evolved stages of high-mass star-forming regions, for instance, the HMC phase (e.g., CH3CN, OCS). The homogeneous list of molecules from various chemical “families” made our search for the best-fit model parameters (T, nH, age) more reliable. The spectra of these molecules were fitted with a Gaussian profile and the corresponding column densities were obtained (see below). For that we excluded all lines that deviate extremely from Gaussian profiles, because this makes the derivation of column densities, based on the assumptions described in Sect. 4.3, unreliable. The reasons for such complex line shapes might be optical depth effects, self-absorption, flux accidentally measured at an 240″ off-position due to the limited distance in the wobbler-switching mode, or peculiar gas kinematics. Furthermore, the spectra of very abundant molecules might be affected by foreground and background emission. Infall or outflow motions also lead to a departure of a line profile from Gaussian shape toward double-peaked profiles. In the source W3(H2O) a masering methanol line caused a high intensity in one channel around the line peak of that line. Other problems in fitting the line profiles were broad line-wings that we were unable to fit with a single Gaussian, but only with a combination of a narrow main component and a broader underlying component. In these cases we fitted the narrow component, which most likely traces the bulk component of the molecular gas. Some of the lines with slightly asymmetric shapes still allowed Gaussian fitting, which we also included in the analysis, with the caution that the derived values are probably underestimated.

The hyperfine line of HCN (10) shows anomalies in the line ratios compared with the predicted theoretical values in a substantial number of sources, which was found by other groups as well (see Cernicharo et al. 1984; Loughnane et al. 2012). In the optical thin LTE case the relative intensities of the three HCN(1-0) hyperfine components are F(01):F(21):F(11) = 1:5:3 (e.g., Wannier et al. 1974). In our spectra we see various deviations from these ratios with either F(01)/F(11) ≧ 1, which cannot be solely explained by the optical thickness of the line. The probable reason for this behavior is scattering of radiation from the core in a moderate-density envelope (Cernicharo et al. 1984). This is possibly enhanced by an inherently complex spatial distribution of HCN through the cloud or clumpiness of the cloud itself, which causes tricky absorption and re-emission of the radiation. In such a case we refrained from using this line in the data analysis. In summary, the main problems in the analysis of the spectral lines were optical depth effects, a poor off-position, foreground and background emission, a complex velocity structure of the source, and anomalies in the HCN hyperfine structure. The assumptions used to derive molecular abundances and column densities are presented in Sect. 4.3, and we discuss the uncertainties in Sect. 4.4.

thumbnail Fig. 1

From left to right, we show sample 16 GHz spectra of the sources IRDC048.6, HMPO18151, HMC029.96, and UCH013.87, representing the four selected evolutionary stages of high-mass star formation. The frequency resolution is ~0.2 MHz and typical rms values are ~0.1 K at 1 mm and ~0.03 K at 3 mm.

Open with DEXTER

Table 2

List of analyzed molecules with transitions, frequencies, energies of the upper level, critical densities, and effective densities calculated for 10 K and 100 Ka.

4. Results

4.1. Chemical characteristics of the four evolutionary stages of high-mass star-forming regions

In Fig. 1 we show sample spectra for each of the four evolutionary stages. We now describe the general characteristics of each stage on the basis of these four examples. In the object IRDC048.6 we only detected simple and common molecules such as isotopologues of CO or HCO+ and tracers of cold and dense matter such as N2H+. The spectra of HMPO18151 already shows a higher number of detected molecules and higher intensities of the observed lines. In this phase the temperature starts to increase, which allows radicals to become mobile on the dust grain surfaces, synthesizing complex (organic) ices and simpler species that may be able to evaporate to the gas phase. This also leads to an overall increase in the intensities of the detected lines (e.g., C2H, H2CO, HCN, HCO+ and isotopologues), such that more species appear on the spectrum. We detected an organic species, methanol (CH3OH), SiO – a molecule that is a tracer of outflows and shocks, and we began to see sulfur-bearing molecules (e.g., CS, SO, H2CS and isotopologues). In the HMC stage the spectrum shows the largest number of lines, including transitions of long and complex molecules such as methyl formate CH3OCHO. Some of the HMC spectra already show recombination lines (e.g. HMC029.96), which is typically only visible in the spectra of UCHii regions. In this case they are generated by a nearby UCHii region (see Beuther et al. 2007b) that is also covered by the single-dish beam. During this late evolutionary phase many complex, photofragile molecules might be destroyed by the intense ionizing radiation of the central star(s). This enhances the amount of ionized medium and thus reduces the amount of the neutral medium compared with HMCs, leading to lower line fluxes and column densities. Correspondingly, the UCHii spectrum shows fewer molecular lines. Apart from additional recombination lines of hydrogen, the UCHii spectra looks qualitatively very similar to the HMPO spectrum.

We selected 15 different molecular species that trace the physical conditions and chemical evolution to analyze and compare them with a chemical model. A list of all the selected lines is given in Table 2. To exclude optically thick lines we used the rarer isotopologues, when available. Specifically, we did not analyze 13CO, HNC, HCO+, or CS. The assumed isotopic ratios are given in Sect. 4.3. The detected lines in each source were fitted with a Gaussian profile to obtain the total line intensities. For sources in which we did not detect a specific transition we used the 3- sigma rms value to estimate an upper limit of the corresponding column density.

4.2. Calculation of the H2 column densities

To derive abundances we first calculated the H2 column densities for all observed regions. Many of our sources are covered by the galactic plane survey ATLASGAL (Schuller et al. 2009), which observed dust emission at 870 μm. Most of the HMPOs were observed at 1.2 mm with the bolometer MAMBO at the IRAM 30m. The data for the remaining sources were taken from the SCUBA Legacy Catalog published by Di Francesco et al. (2008) which are bolometer maps observed at 850 μm with the JCMT. In Table 1 we list which continuum data were taken for each of the sources.

The coordinates for the spectral observations were centered on the highest intensity peak in the bolometer map. With the observed value of the peak intensity we derived H2 column densities following Eq. (1) taken from Schuller et al. (2009). The typical temperatures for the observed sources were assumed. Average IRDC temperatures are around 15 K based on several NH3 IRDC surveys (e.g., Sridharan et al. 2005; Pillai et al. 2006; Chira et al. 2013). Typical HMPO temperatures were estimated from SED fits and are ~50 K. HMCs are usually warmer and we chose T = 100 K as the average temperature. The same temperature of T = 100 K is assumed for the UCHii stage. Finally, the dust opacities were interpolated from Ossenkopf & Henning (1994), assuming grains with thin ice mantles, gas densities of n = 105 cm-3, and a gas-to-dust mass ratio R = 100. The actually used dust opacities are κ850 μm = 1.48, κ870 μm = 1.42, κ1.2 mm = 0.97. Using these assumptions and assuming that the emission is optically thin and at LTE, the H2 column density is calculated as (1)In this calculations we also assumed that the dust and the gas are collisionally coupled and thus have the same temperature. The observations with the IRAM 30 m at 1 mm have a HPBW of 11″ and at 3 mm a HPBW of 29″. The data obtained with MAMBO have a resolution of 11″, the ATLASGAL maps 19.2″ and the SCUBA data 22.9″. For the molecules observed at 3 mm we smoothed the maps with a Gaussian kernel to obtain the same resolution and derived the corresponding H2 column density. For the molecules observed at 1 mm we used the original dust continuum maps to calculate the H2 column density. Thus we have two different H2 column densities measured for each source on two different spatial scales. While the H2 column densities for the molecules at 3 mm are beam matching, this is not the case for all sources for the molecules at 1 mm. We slightly underpredict the H2 column densities in the cases of the 1 mm molecular line data, where we only have ATLASGAL or SCUBA data because of the larger beam. To facilitate the modeling, we calculated for each angular resolution the corresponding physical radius, using estimates of the distances to the sources. Thus, we avoided effects from the heterogeneous nature in angular resolution of the continuum and spectroscopic data. For sources with a distance ambiguity we used the lowest value. Three of the sources show different velocity components in some transitions, resulting in different estimates of the kinematic distances. In these cases we chose the kinematic distance corresponding to the main velocity component that was detected in all transitions.

Another uncertainty in the comparison of the continuum data and the molecular line data comes from the different techniques of measuring. The continuum data are observed with bolometer arrays and suffer for some degree from spatial filtering of the large-scale structures (see Motte & André 2001; Schuller et al. 2009). In contrast, the spectroscopic data are observed with heterodyne receivers and have no filtering problems. This is a general uncertainty in all comparisons between bolometer- and heterodyne-receiver-based measurements and leads to an underestimation of the true flux of the bolometer-based measurements. Hence, in our case the derived abundances are slightly overestimated. This effect is hard to quantify, but is not expected to play an important role in our case. It is reasonable to assume that the large-scale structure, which is filtered out, consists of diffuse gas with low density. A possible way to quantify this effect of missing large-scale structures in a first approximation is to take the canonical CO abundance of 10-4 and compare it with our derived values, which are about a factor 24 lower. Thus, the estimated H2 column densities are even higher than the canonical values and the spatial filtering seems to have only a minor impact.

Table 3

Observed median abundances and standard deviation in a(x) = a × 10x for the IRDCs for T = 15 K and T = 20.9 K. x(e) is the ionization fraction.

Table 4

Observed median abundances and standard deviation in a(x) = a × 10x for the HMPOs for T = 50 K and T = 29.5 K. x(e) is the ionization fraction.

Table 5

Observed median abundances and standard deviation in a(x) = a × 10x for the HMCs for T = 100 K and T = 40.2 K. x(e) is the ionization fraction.

Table 6

Observed median abundances and standard deviation in a(x) = a × 10x for the UCHiis for T = 100 K and T = 36.0 K. x(e) is the ionization fraction.

4.3. Calculation of the molecular column densities

To calculate column densities for the different molecules in Table 2 we made several simplifying assumptions. The uncertainties introduced by these assumptions are discussed in Sect. 4.4. 1). Local thermodynamic equilibrium (LTE), since the densities in the observed high-mass star-forming regions are on the order of 105 cm-3 and higher, which exceeds the critical densities for many of the observed transitions; 2) a uniform gas kinetic temperature for all molecules observed in the same evolutionary stage. In the first iteration of the comparison with the model we assumed T = 15 K for IRDCs, T = 50 K for HMPOs, T = 100 K for HMCs and T = 100 K for UCHiis, respectively (iteration 0). In a second step we changed the assumed temperatures based on the outcome of the model (iteration 1). Due to the large beam sizes used to observe the high-mass star-forming regions, we took lower average temperature values for the later stages than the values observed with high-density tracers such as methyl cyanide; 3) the dust and gas temperatures are equal; 4) line emission was assumed to be optically thin. For the species for which a rarer isotopologue line is also detected, we used this line to derive the column density (in particular, HN13C, H13CO+, 13CS, C33S, and C18O). For N2H+ and HCN we fitted the hyperfine line structure and constrained the optical depth; 5) derived column densities were beam averaged; 6) the medium is spatially homogeneous; 7) the Rayleigh-Jeans approximation of the Planck function was used; 8) no isotopic fractionation. Because we only provide single-dish measurements with limited spatial resolution at single pointings, only mean values smoothed over the beam were derived.

We calculated the molecular column densities following the equations in Tielens (2005): (2)where the line frequency νul is in GHz, the integrated intensity is in K km s-1, and the Einstein coefficient Aul is in s-1. Then the total column density can be calculated, (3)where Q is the partition function, gu the statistical weight of the upper level, Eu the upper state energy, and k the Boltzmann constant. Tex is the excitation temperature and assumed to be equal to Tkin. Finally, we converted the derived column densities for the rarer isotopologues given in Table 2 into the column densities of their main isotopologues. We hereby assumed for all sources the same relative isotopic ratios representative of the Sun and local ISM: 12C/13C = 89, 16O/18O = 499, and 32S/33S = 127 (Lodders 2003). The molecular column densities were then divided by the H2 column densities and the abundances were derived. The derived abundances for the molecules with transitions at 1 mm are upper limits because of the differences in resolution between ATLASGAL/SCUBA and the 1 mm molecular data. This does not affect the modeling, because we directly converted the beam sizes into physical sizes. The resulting median abundances including all detections and upper limits for each subsequent evolutionary phase are given in Tables 36. In addition to the molecules listed in Table 2, in some regions we detected also HC3N, H2CS, CH3CHO, and other species. However, analyzing these is not a main focus of the current paper. The more detailed investigation of these species will be presented in future publications.

4.4. Uncertainties in the observed values

Our estimates of the column densities are based on simplifying basic assumptions, such as optically thin line limit, uniform excitation temperature for all observed molecules, and all sources being at a particular evolutionary stage, uncertainties in the derived H2 column densities, which all need to be taken into account in the data analysis and further modeling.

We assumed a fixed dust opacity for all sources. In addition, we assumed a single excitation temperature for all sources and molecules in one evolutionary sample that is equal to the assumed kinetic temperature in this evolutionary stage. In this way, we did not account for source peculiarities or a possible misclassification of the stage or an overlap of several stages in a single source. Furthermore, the kinetic temperature of single molecules might deviate from the real excitation temperature in case of subthermally excited lines. This might especially be the case for the HMCs and UCHii regions, which have the highest kinetic temperatures. Comparing the calculated abundances derived with the high excitation temperature Tex = 100 K with those derived with the lower excitation temperature of Tex = 50 K, we found that the abundances decrease by factors of 3 − 10 in the case with the low Tex. HNCO, CH3CN, CH3OH, and CH3OCHO are the most affected molecules.

We made the same exercise for the IRDCs and compared the abundances calculated with Tex = 15 K with those calculated with 10 K. The difference in the abundances reaches at most a factor of 3, with the exception of OCS, which shows a large difference of a factor of 15. We attributed this deviation to a non-LTE effect because a high J rotational line of OCS was employed, which might be subthermally excited even at typical densities of the high-mass star-forming regions. A more detailed discussion of the effect of uncertainties on the excitation temperature is given in Sect. 5.3. Furthermore, assuming optically thin lines implies that the calculated column densities may also be just lower limits.

Moreover, uncertainties in the conversion of the dust into gas densities constitute a factor of about 1.5 compared with Draine (2011). The isotopic ratios vary between different environments and the uncertainties in the adopted isotopic ratios are a factor of 2 − 4 (e.g., Wilson & Rood 1994). Glassgold et al. (1985) found that the 12C/13C ratio in PDRs might deviate even by more than a factor of 10 from the solar value. More recently, Röllig & Ossenkopf (2013) suggested strong fractionation effects in C-chemistry.

In combination, these systematic errors result in an overall uncertainty of about one order of magnitude in the derived molecular abundances. This level of uncertainty is typical for astrochemical studies. It is present in the assumptions made to analyze observational data as well as in the intrinsic uncertainties of chemical models (e.g., Wakelam et al. 2010; Vasyunina et al. 2011, 2012; Albertsson et al. 2013). Improvements of this situation are an area of active study.

In this work, these uncertainties were taken into account during the data fitting, as explained in Sect. 5.1.4. To avoid suffering strongly from possible misclassifications and peculiarities of single sources, we treated the sample in a statistical way and were mainly interested in the median characteristics of a subsample. That might lead to larger spreads within one subsample, but makes the analysis and comparison with the model more robust and reliable.

thumbnail Fig. 2

Relative detection fractions for each analyzed molecular transition ordered from left to right: IRDCs without associated point sources at wavelengths below 70 μm, IRDCs with associated point sources, HMPOs, HMCs, and UCHiis. The color notation for the different source types is explained in the upper right corner.

Open with DEXTER

4.5. Chemical evolution

In Sects. 4.2 and 4.3 we calculated the column densities of H2 and other molecular species. From these values we derived averaged abundances for each source. In the next step we combined all these abundances from all sources at a specific evolutionary stage and derived a characteristic chemical “portrait” for each phase. The results are given in Tables 36 together with the detection fractions. The detection fractions are plotted in Fig. 2 and the abundances are plotted in Fig. 3 for the IRDCs, HMPOs, HMCs and UCHii regions. The individual column densities for each stage are presented in Figs. A.1A.4.

All abundances in all four evolutionary stages have values ranging between 10-11 and 10-6. The spread of column densities within a particular stage is about one order of magnitude or lower, except for H2CO in the IRDCs and N2H+ in the HMCs with spreads of almost three orders of magnitude. Within the two subsamples of IRDCs (with and without associated point sources at wavelengths below 70 μm) the difference in median abundances does not exceed a factor of 2. But in terms of detected molecules we found a stronger difference. The molecules 13CS, CH3CN, and OCS are only detected in the more evolved IRDCs. The detection fractions of SiO, SO, and CH3OH are at least 20% lower toward IRDCs without a point source than in sources with an embedded point source. This is a sign that the latter group is already at a slightly more advanced chemical evolutionary stage.

In general, the detected abundances over all evolutionary stages tend to increase with rising temperature during the evolution. To identify the chemical traces of the evolution we compared the abundances within the 2575% range of the median value between the different evolutionary phases. The abundances in the IRDC stage are around 10-10 and 10-9, with OCS, C2H, and the CO isotopologue having abundances higher by one to two orders of magnitudes. Compared with the HMPOs, all abundances in the IRDCs are lower by about a half to one order of magnitude. Exceptions are the organic species CH3OH, CH3CN, and H2CO, which have abundances similar to those in the HMPOs. Moreover, OCS is remarkably much more abundant, but only detected in IRDC028.2. The overall detection rate in the IRDCs is lower than that of the HMPOs. Next, in the HMPOs the abundances have values mainly in the range of 10-10 to 10-8. They are even higher in the HMCs. Only N2H+ has similar abundances across these stages. The abundances of H13CO+ and HN13C, while still higher in the HMC phase, are closer to the respective HMPO values than the values of other molecules. The abundances of most of the molecules in the HMCs have values between 10-9 and 10-7. This is the only phase where we detected the most complex molecule analyzed in our sample, methyl formate CH3OCHO. However, methyl formate was detected in IRDC028.2 by Vasyunina et al. (2014), which means that this non-detection might be a combination of lower sensitivity and a weaker transition.

From the HMCs to the UCHii stage, abundances of CH3OH, CH3CN, OCS, C33S, H2CO, SO, and SiO become lower. The remaining molecules have similar or slightly higher abundances in the UCHii stage. However, the overall detection rate of the different molecules for the HMCs is much higher than for the UCHiis. Complex or heavy molecules such as CH3OCHO, CH3OH, SiO, HNCO, and OCS are not detected in 30% of the UCHiis and are almost always found in the HMCs.

thumbnail Fig. 3

From top to bottom abundances (with respect to the total H2) of the analyzed molecules in the IRDC, HMPO, HMC and UCHii sample. The red line shows the median, the green cross is the mean, the bar indicates the inner 25%75% range around the median and the whiskers mark the total range of all detections. The black arrows indicate the lowest upper limit of all calculated upper limits for this particular molecule and stage.

Open with DEXTER

4.6. Ionization degree

In our spectral setup we observed two molecular ions, H13CO+ and N2H+. Together with the observed H2 densities this enables us to calculate lower limits of the electron fraction following Caselli et al. (2002). This method is based on calculations from Caselli et al. (1998) using a simple chemical model and is only valid for homogeneous clouds. Since the deuterated counterparts are not covered by our setup we were only able to derive the electron fraction based on these two molecular species. The median electron fractions are given in Tables 36. The median value increases with evolutionary phase from 5 × 10-9 in the IRDCs to ~10-8 − 10-7 in the UCHii regions (depending on the assumed excitation temperature). This behavior agrees with the scenario of an increasing ionization radiation during the evolution of the central source that finally reaches its maximum in the UCHii phase. The derived values are also consistent with a lower limit from Caselli et al. of  10-9 and the work by Miettinen et al. (2011), who found lower limits of ~10-8 − 10-7 in massive clumps associated with IRDCs.

5. Discussion

In the previous section we described the observational results about the chemical composition of the different evolutionary phases. In the next step we model the chemistry over the full evolutionary range. A detailed background of the presumed evolutionary sequence is outlined in the introduction in Sect. 1. Furthermore, combining the information from observations and model, we give a comprehensive description of this evolutionary sequence in Sect. 5.5.

5.1. Model

In this section we describe our iterative chemical fitting model “MUSCLE” (“MUlti Stage CLoud codE”), which we used to constrain mean physical properties and chemical ages during the evolution from the earliest IRDC phase till the late UCHii phase of the high-mass star formation.

5.1.1. Physical model

We modeled each evolutionary stage in a 1D approximation. Each environment is spherically symmetric, homogeneous, and has a fixed outer radius of rout = 0.5 pc. This radius is close to a typical value for dense parts of high-mass star-forming regions and also represents the largest 29″ IRAM beam size used in our single-dish observations as derived for the mean distance of the sample. The radial density and temperature profiles are modeled by modified power laws, (4)and (5)respectively. Here, rin is the inner radius below which temperature and density are assumed to be constants and p and q are the corresponding radial profiles of density and temperature. The parameters rin, p, ρin and Tin are the fitted physical quantities. rin was varied between 10-4 − 10-2 of rout. The gas and dust were assumed to be in thermal equilibrium, therefore we did not calculate dust and gas temperatures separately. Equation (4) was numerically integrated from r = 0 till r = rout to derive the total modeled H2 column density of a cloud. For that, a grid of 30 − 40 radial cells was found to provide sufficient accuracy, and hence was adopted in large-scale chemical modeling.

We assumed that at each stage a cloud is furthermore embedded in large-scale, low-density matter that shields it from the interstellar FUV radiation by visual extinction of 10 mag. For UCHii regions we assumed that the size of a forming Strömgren zone is still small, 0.1 pc, such that it can be neglected for calculations of bulk beam-averaged quantities such as molecular column densities.

5.1.2. Chemical model

The adopted time-dependent gas-grain chemical model based on the “ALCHEMIC” code is fully described in Semenov et al. (2010). A brief summary of the updates is provided below. The chemical rate file is based on the 2007 realization of the OSU network2. The network is supplied with a set of ~1000 reactions with high-temperature barriers from Harada et al. (2010) and Harada et al. (2012). The recent updates as of September 2012 to the reaction rates are implemented (e.g., from KIDA database3), see Albertsson et al. (2013). We considered cosmic-ray particles (CRP) and CRP-induced FUV radiation as the only external ionizing sources. The CRP ionization rate ζCR = 5 × 10-17 s-1 was used. The UV dissociation and ionization photorates from van Dishoeck et al. (2006)4 were adopted, assuming the case corresponding to the spectral shape of the interstellar FUV radiation field.

In addition to pure gas-phase chemical processes, the chemical network includes gas-grain interactions. We assumed that molecules stick to grain surfaces at low temperatures with 100% probability. We did not allow H2 to stick to grains because the binding energy of H2 to pure H2 mantle is low, ~100 K (Lee 1972). The ices are released back to the gas phase by thermal, CRP-, and CRP-induced UV-photodesorption. The grain re-charging was modeled by dissociative recombination and radiative neutralization of ions on grains, electron sticking to grains. Chemisorption of surface molecules was not considered. We used the UV photodesorption yield for ices of 10-3 (e.g., Öberg et al. 2009a,b).

The synthesis of complex molecules was included using a set of surface reactions (together with desorption energies) and photodissociation reactions of ices from Garrod & Herbst (2006) (see also Semenov & Wiebe 2011). We assumed that each 0.1 μm spherical olivine grain provides 1.88 × 106 surface sites for accreting gaseous species. The surface recombination solely through the Langmuir-Hinshelwood formation mechanism is considered. Upon a surface recombination, the reaction products are assumed to remain on the grain as the grain lattice would absorb all energy released during the recombination. Following interpretations of experimental results on the formation of molecular hydrogen on dust grains (Katz et al. 1999), we employed the standard rate equation approach to the surface chemistry without H and H2 tunneling through the potential walls of the surface sites. Overall, our chemical network consists of 656 species made of 13 elements, and 7907 reactions. No chemical parameter was varied during the iterative fitting of the observational data.

Our ALCHEMIC Fortran code (Semenov et al. 2010) is based on the double-precision variable-coefficient ordinary differential equation solver with the preconditioned Krylov (DVODPK) method GMRES for the solution of linear systems5. The approximate Jacobi matrix is generated automatically from the supplied chemical network. For astrochemical models dominated by hydrogen reactions the Jacobi matrix is sparse, with ≲ 1% of non-zero elements. The corresponding linearized system of algebraic equations is solved using a high-performance sparse unsymmetric MA48 solver from the Harwell Mathematical Software Library6.

5.1.3. Initial abundances

The initial abundances very close to the low-metals set of Lee et al. (1998) were used to model the chemistry at the IRDC evolutionary stage (see Table 7). The only difference was the initially atomic abundance of Si (3 × 10-9 with respect to H) and S (8 × 10-7 with respect to H), which we had to use to achieve a better fit to the observational data. Other initial abundance sets with different amounts of oxygen and sulfur, which have been proposed in studies of chemistry in molecular clouds and protoplanetary disks, did not produce as good a fit to our IRDC data (see e.g., Wakelam & Herbst 2008; Graedel et al. 1982; Jenkins 2009; Hincelin et al. 2011; Dutrey et al. 2011). For modeling the chemistry of the evolutionary stages after the IRDC stage, the final chemical abundances from the best-fit model of the previous evolutionary stage were used as initial abundances (i.e., to model HMPOs we used the molecular abundances of the best-fit IRDC model at the best-fit chemical age). This allowed us to reproduce an approximately steady warming of the matter during the evolution of the high-mass star-forming clouds.

Table 7

Initial atomic and molecular abundances.

5.1.4. Iterative fitting of the data

Using the physical and chemical models described above, we iteratively fitted the mean observed molecular column densities for each individual evolutionary stage by varying the inner radius rin, the density at the inner radius ρin, the temperature at the inner radius Tin and the density slope p. First, for each 1D cloud model with a grid of 30 cells we derived appropriate ranges of rin, the ρin, and the density profile ρ by calculating beam-averaged H2 column densities such that they are within a factor of 3 to the observationally derived values. To limit the computational time and storage space required for our large-scale modeling, we restricted the density profile parameter p to lie between the values of 1.5 and 2.0, as seen in a variety of observations by Beuther et al. (2002a), Mueller et al. (2002) and Hatchell & van der Tak (2003). The inner radius was varied between 5 × 10-5 − 5 × 10-2 pc. Another parameter whose value was iteratively varied was the temperature at the inner radius, Tin (between 15 − 27 K for IRDCs, 50 − 250 K for HMPOs, 150 − 500 K for HMCs and UCHii regions). All other physical parameters were fixed during iterations, including the temperature profile q, which we set prior to the modeling. In the first stage we assumed an isothermal sphere with q = 0 and in the later stages a temperature profile with fixed q = 0.4 as a standard value. Simulations of massive star formation by van der Tak et al. (2000) found that the temperature profiles have steeper indices in the inner regions, but flatten outside the inner 2000 AU and run asymptotically into a distribution with q = 0.4. Typically, by varying the parameters rin, ρin, Tin and p for each stage, one model realization consists of 5003 000 cloud models.

Using this large set of cloud models and our chemical model described above, we calculated corresponding time-dependent chemical structures over 1 Myr (with 99 logarithmically taken time steps between 104 and 106 years). After that, we calculated the corresponding beam-averaged column densities and the goodness of the fit for each cloud model, observed molecule, and each time step by using a confidence criterion from Eq. (3) in Garrod et al. (2007). As the standard deviation σ for all observed molecules we assumed an order of magnitude variation between the modeled and the observed column densities. For molecules for which only the upper limits of the observed column densities were available in more than 50% of the sources (due to non-detections), we set the goodness-of-fit to be 0 when the modeled column density was lower than the upper limit multiplied by 10. Then, the confidence criteria for all molecules at a given time step and for a given cloud model were summed and divided by the number of observed molecules, and a total goodness-of-fit was obtained. Finally, a code searched for the lowest value that represented the best-fit physical model and the best-fit chemical age. Because we modeled averaged data, the chemical age is also an average for each of the evolutionary stages. These ages can be considered as a typical mean time for an object with a certain chemical composition to evolve to the stage with the subsequent characteristic chemical composition. A typical execution time of the MUSCLE model for one evolutionary stage is between 2 and 10 h on a server with 24 Intel Xeon 3.2 GHz CPUs.

Table 8

Parameters of the best-fit IRDC model.

5.2. Radial distributions of molecules in the model

thumbnail Fig. 4

Snapshot of the radial distributions of CO, HCO+, N2H+, C2H, H2CO, CH3OH, CH3CN, and SiO at the time of the best-fit model of the HMC phase. The time is given relative to the beginning of the HMC phase.

Open with DEXTER

In Figure 4 we show the radial distribution of several molecules for the best-fit model of the HMC stage. Simpler molecules like N2H+ and HCO+ , which are produced in the gas-phase, have a flatter distribution or even increase with radius, whereas more complex molecules like CH3OH and CH3CN are predominantly abundant in the warm inner region, where ices are able to thermally desorb from dust grain surfaces. This leads to the conclusion that to properly account for these effects it is also necessary to assume different beam-filling factors for different molecules. Interferometric observations with high spatial resolutions are needed to address this problem in detail.

5.3. Comparison of observations and the best-fit models

The comparison of the data with the best-fit model described in Sect. 5.1.4 for our samples of IRDCs, HMPOs, HMCs, and UCHiis, shows good overall agreement between the observed and modeled column densities. The best-fit parameters are shown in Tables 811 and their corresponding column densities in Tables 1215. Furthermore, we list the corresponding radial temperature and density profiles as well as the chemical age of the individual stages derived with the time-dependent chemical code. Figure 5 shows the observed and modeled abundances as a function of evolved time. In this plot, the modeled abundances between the different stages seem to have discontinuities. This is because the best-fit models of various stages do not have similar physical structures. The chemistry reacts quickly to these changes in the physical structure at the beginning of a new stage. The abundance at each subsequent stage is plotted 100 years after the starting time of the new stage and thus the abundance occurs as a discontinuity.

Table 9

Parameters of the best-fit HMPO model.

Table 10

Parameters of the best-fit HMC model.

Table 11

Parameters of the best-fit UCHii model.

Table 12

Median column densities in a(x) = a × 10x for observations and best-fit IRDC model.

Table 13

Median column densities in a(x) = a × 10x for observations and best-fit HMPO model.

Table 14

Median column densities in a(x) = a × 10x for observations and best-fit HMC model.

Table 15

Median column densities in a(x) = a × 10x for observations and best-fit UCHii model.

In the IRDC stage the best-fit model was able to successfully explain 14 out of 14 different species. The best-fit age is ~11 000 years (with an uncertainty of a factor of 23). The chemical age of a given stage is the time the model needed to evolve from the starting conditions to the best-fit model of this stage. For the later stages, the chemical composition of the best-fit model of the previous stage was set as the starting conditions for the next one.

In the second stage of an HMPO the best-fit model was able to successfully explain 12 out of 14 different species. The two species that are not explained are H2CO and CH3OH. The model strongly overproduces H2CO by a factor of more than 100. Methanol is underproduced, but only slightly lower than the assumed criteria of one order of magnitude. In general, problems with fitting the data are due to environmental effects that are not taken into account in the model, such as shocks and UV-penetration. CH3OH and SiO are known shock tracers and thus more likely to be not explained by the model. A not yet fully understood surface chemistry is another source of uncertainty in the model. H2CO is partly produced in the gas phase and on the grains. The two different formation routes are strongly time dependent, and so far there is no consensus in the community about how to implement this dependency or the implementation of different desorption mechanisms that allow molecules to return to the gas phase (see e.g., Garrod et al. 2007; Vasyunin & Herbst 2013b). In this case the timescale seems to be too short to transform H2CO into CH3OH, which can partly explain the overproduction and underproduction, respectively. The best-fit age of the HMPO phase is ~60 000 years (with an uncertainty of a factor of 23).

In the third stage of an HMC the best-fit model was able to successfully explain 11 out of 14 different species. The three species that differ are HNC, CS, and C2H. While HNC and CS are only slightly underproduced compared with the agreement criteria, C2H is underproduced by almost more than a factor of 50. In the model the HMC is still surrounded by an external shell with a visual extinction of 10 mag with no external UV-field, which leads to lower C2H abundances. In reality, the cloud might be highly fractal and fragmented such that UV-photons can penetrate deeper and produce more C2H. The best-fit age of the HMC phase is ~42 000 years (with an uncertainty of a factor of 23).

In the fourth stage of an UCHii the best-fit model was able to successfully explain 11 out of 14 different species. The three species that are not reproduced are HNC, SO and C2H. A caveat for our models of UCHii regions is that it has an internal UV-source that is not properly taken into account in the model. Furthermore, since this stage depends on all best-fit models of the previous stages, it is the most doubtful one. The successful fitting of 80% of the molecules is very remarkable already. The best-fit age of the UCHii regions is ~11 500 years (with an uncertainty of a factor of 23).

This is the first time that the chemical evolution was modeled and reasonably well fitted over the entire span of the early evolution of high-mass star formation regions. The mean temperatures of the best-fit model temperature profile for the last three stages (see Tables 811) are considerably lower, by a factor of 23, than the excitation temperatures used to calculate the observed column densities. This difference is much smaller in the IRDC stage.

thumbnail Fig. 5

Relative abundances to H2 plotted for all four stages. The modeled values are shown as a black solid line, the observed values are depicted by a blue dashed line. The error bars are indicated by the vertical marks.

Open with DEXTER

In a second step we took the mean temperatures from each best-fit model as excitation temperatures to recalculate the observed column densities. The difference in abundances for the two different excitation temperatures for the IRDCs is at most a factor of 2 for all molecules, except for OCS, which differs by a factor of 3. The abundances in the HMPOs change at most by a factor of 3. The discrepancy between the temperatures originally used and the mean temperatures derived in the best-fit for the HMCs and UCHii regions is larger than for the IRDCs and HMPOs. Thus, the change in the abundances is higher by about a factor of 6 for the HMCs (only CH3OH and CH3OCHO differ by factors of 7 and 10, respectively) and 69 for the UCHiis (only HNCO and CH3OCHO differ by a factor of 13). These values are shown in Tables 36 as well.

In this iterative step we sought again to find best-fit solutions for the newly derived values with more consistent temperatures. To keep the model physically meaningful, the temperature profile of the last three stages was set as a power-law with a constant temperature of the IRDC best-fit value in the outer parts. This resulted in a stable solution for the IRDCs. But for the later stages the model tends to produce a small and dense inner hot part and a large isothermal outer part. A decrease in the modeled best-fit temperature leads to higher observationally derived H2 and lower observationally derived molecular column densities. Fitting the new higher H2 and lower molecular column densities lowers the mean temperatures even more and leads to a non-converging feedback. However, temperature measurements of tracers of the inner dense part such as CH3CN show a hot inner component with T > 50 K.

The discrepancy in temperatures between model and observations is caused by uncertainties and simplifications in both. On the observational side we have single-pointing single-dish spectra for which we assumed a single excitation temperature. However, the observed sources are most likely inhomogeneous and have structures of differing environments. We partly took this into account in the model by assuming a 1D physical structure with a temperature and density profile. Another reason for this mismatch could be subthermal excitation. To avoid the effect of non-converging feedback and to account for all the caveats, we stopped the modeling process and report the results of the first iteration of the best-fit model.

Based on the probable explanations for the mismatch between the temperatures used for the observations and derived with the model, improvements on both sides might help to solve the discrepancy. Detailed gas and dust temperature maps and dust continuum radiative transfer calculations might solve the degeneracy in constraining the temperature. Interferometric observations with higher resolutions would help to decrease uncertainties in the beam filling fraction. Moreover, extending the observed species toward deuterated molecules might help to better constrain temperatures since deuterium chemistry is sensitive to the previous thermal history of the environment (see Caselli & Ceccarelli 2012). In addition, the model might be improved by including radiative transfer, a gradual warm-up phase (see Garrod et al. 2008; Vasyunin & Herbst 2013a), and in general more physics like the proper treatment of shocks, outflows and UV radiation transport. Furthermore, at least 2D models would be necessary to be able to explain observations with higher resolutions.

5.4. Comparison of observed and best-fit molecular ratios

In addition to the information from column densities of single molecules, the ratios between different molecules provide another interesting tool for investigating the chemical evolution. In Fig. 6 the observed and best-fit ratios of several molecular pairs are shown. Due to the limited sensitivity in the observations the overall dynamic range of observed ratios is smaller than in the model. In general, the ratios of the observed median column densities show no obvious trends with evolution that cannot be explained within the uncertainties. Because of the poorly explained molecular column densities in the best-fit model for a few species, such as H2CO or C2H, and the agreement criteria between model and observations, some of the ratios deviate strongly, depending on the specific molecule and evolutionary phase.

5.4.1. HCN/HNC

Previous observations showed that the ratio of HCN/HNC depends on the temperature of the medium, for example, observations of Orion-KL show a decreasing ratio from the warm inner core toward the cold outer edge from 80 to 5 (Goldsmith et al. 1986; Schilke et al. 1992). The best-fit model follows this trend and increases until it reaches the highest value of ~50 in the beginning of the HMC phase, whereas our observations do not show significant differences between the four stages with a ratio between 0.3 − 0.6. The total value of the observed ratio <1 might be underpredicted due to the assumptions of the isotopic ratio in the conversion from HN13C to HNC. Furthermore, due to self absorption, high optical depth, and other effects discussed in Sect. 3.2, we refrained from using the observed HCN lines in the analysis in many cases. In these cases, the sources more likely have higher column densities. That might bias the derived column densities toward lower optical depths and therefore lower column densities.

5.4.2. HCO+/N2H+

The ratio of the two molecules N2H+ and HCO+ is a candidate for a chemical clock. There is evidence that their abundances are anticorrelated in chemical models of low-mass star-forming regions (e.g., Jørgensen et al. 2004) and massive star formation (Barnes et al. 2013), since both molecules are sensitive to the amount of CO. According to the models, N2H+ is more abundant during the early phases when CO is frozen out and is destroyed after CO returns into the gas phase. This behavior of the anticorrelation between CO and N2H+ is seen by Bergin et al. (2002) toward the low-mass star-forming region B68 and in the ratio of N2H+ and CO of this work. The abundance of HCO+ is directly dependent on the CO abundance and the ionization degree and thus is anticorrelated with N2H+. Our observations reproduce this trend of an increasing ratio with evolutionary stage. However, in the model the ratio is almost constant within the given errors. A possible explanation are the uncertainties in the observed values, which are larger than the total change of that particular ratio. Nevertheless, observations and model agree within the uncertainties and show a total spread of the ratio between 10 and 100.

5.4.3. CH3OH

In our model, the ratios of the more complex molecule methanol with simple molecules such as CO is almost constant at 10-6 during the IRDC and HMPO phase and jumps at the beginning of the hot core phase to 10-3. During the IRDC and HMPO phases CH3OH is steadily produced on and in ice mantles of dust grains. One of our model restrictions is that the methanol ice can only be desorbed to the gas by cosmic-ray particles (CRP) and CRP-UV heating. In the HMC phase the entire solid reservoir of dust mantles is released into the gas phase, leading to a sudden increase of the gas-phase concentration of many complex molecules. Their abundance may decrease with time after that event if they are actively destroyed or converted to other species in the gas, because surface chemistry becomes insignificant at 100 − 150 K. Afterwards it steadily declines until it reaches 10-5 at the end of the evolution. This behavior is imprinted by the abundance variability of CH3OH and is thus seen in most of its ratios. It agrees with the observational ratio for the last three stages, but underestimates the IRDC stage.

5.4.4. S-bearing molecules

The ratio of the S-bearing molecules CS and SO is poorly described by the model and represents a well-known problem with sulfur chemistry in modern astrochemical databases (Druard & Wakelam 2012; Herpin et al. 2009; Loison et al. 2012; Dutrey et al. 2011). While in the observations for the IRDC stage both molecules have only an upper limit, the later stages show an increasing trend of the ratio between 10 and 100. In contrast, the model predicts a strongly decreasing CS/SO ratio over about three orders of magnitude. The ratio of SO with the early type molecule N2H+ is also poorly described by the model. Only the HMPO and HMC phase agree slightly, with a decreasing ratio with time. Since the model for CS and SO alone does not agree very well for both, it is natural that the ratio does not agree either.

For completeness and to outline the strengths and caveats of our modeling approach, Fig. 6 presents a few more ratios that we do not discuss in detail here.

thumbnail Fig. 6

Column density ratios for 15 different combinations of molecules plotted for all four stages. The modeled values are shown by the black solid line, the observed values are depicted by the blue dashed line. The error bars are indicated by the vertical marks.

Open with DEXTER

thumbnail Fig. 7

Averaged relative abundances to H2 for electrons, HCO+, and N2H+ plotted for the all four stages. The observed values are shown as solid lines, the modeled values are depicted as dashed lines. The electron abundances are plotted in red, N2H+ abundances in black, and HCO+ abundances in blue. The error bars for both the model and the observational data are a factor of about 3 (not shown).

Open with DEXTER

5.5. Chemical evolutionary sequence – synthesis of observations and models

With observational and modeled data at hand we can derive a general picture of the chemical evolution in high-mass star-forming regions. The starting condition at the year 0 in our model is a cold dense molecular cloud, which is still close to isothermal with the initial chemical composition given in Table 7. Any evolution prior to the IRDC stage to these initial conditions is not considered in this work. In our IRDC model the temperature remains constant, and the chemical evolution proceeds under these conditions for about 11 000 years. At this time the best overall agreement between the observed and the calculated molecular column densities is reached. The modeled abundance of the best-fit model at the best-fit age are used as input abundances for the next HMPO stage of the evolution. From that time the HMPO phase begins in our model. An internal source heats the surrounding medium, such that the temperature increases toward the inner region. This is implemented in the modeling by a temperature profile q = 0.4. Increased temperature results in a higher mobility of the surface molecules and a more efficient desorption of ices into the gas phase. Consequently, the molecular abundances increase. Our best-fit model matches the observed values for the HMPOs after 60 000 years (with a total age of 71 000 years from time zero). As temperatures and densities continue to rise due to formation of a central protostar(s), the object enters the HMC phase. The subsequent increase in temperature releases most of the molecules locked along with water ice in the grain mantles into the gas phase. This increases gas-phase abundances of already abundant molecules even further and observed transitions become stronger. The detection rate of molecular lines increases as well, including those from complex and heavy molecules. According to our best-fit model, the HMC phase lasts for about 42 000 years, leading to a total age of 113 000 years. Finally, the last considered stage of high-mass star-forming regions begins. For the UCHii phase the temperature and density structures are re-adjusted again. The best-fit age derived for this last stage is 11 500 years. Thus, we predict that the total evolution from the formation of an IRDC until the formation of an UCHii region takes about 124 500 years (with an uncertainty of a factor of 23). It appears to be consistent with the estimates obtained from modeling of the formation of high-mass stars (e.g., McKee & Tan 2003).

5.6. Observed ionization degree compared with model data

In Sect. 4.6 we calculated lower limits for the ionization degree in the four evolutionary stages of high-mass star formation. In addition, the best-fit model provides all necessary information to derive the exact modeled ionization degrees to compare. The modeled abundances of the main constituents, namely, electrons, N2H+, and HCO+ and the observed abundances for N2H+ and HCO+, are shown in Fig. 7. While the observed lower limit increases with evolutionary phase, the ionization degree in the best-fit model remains almost constant at ~10-7 throughout the whole evolution. In comparison with the observed data it is off by about one order of magnitude in the first two stages and agrees for the last two stages. This leads to the conclusion that a lower limit for the ionization degree derived from HCO+ and N2H+ is a good first approximation. However, it also indicates that in particular in the early evolutionary phases, other ions contribute more strongly to the total ionization degree. These are deuterated ions in the dense cores, which we will observe in a follow-up study, but the singly ionized carbon is also important during these stages. For example, Miettinen et al. (2011) took deuterated ions into account as well and found slightly higher lower limits for x(e)~ 10-8 − 10-7, which are comparable with the results of our model. In addition to this, more sophisticated methods are needed to estimate the ionization degree with an uncertainty better than one order of magnitude. Here, any deviation from the assumption of a homogeneous medium used in our model (e.g., clumpiness on various scales) will tend to increase the modeled ionization fraction for the later evolutionary stages when some high- energy photons might already be produced by forming protostar(s).

5.7. Comparison with literature

Our comprehensive study of the chemical evolution in the early phases of massive star formation can be set in context to previous studies focusing on distinct evolutionary phases. The comparison of our work with the study of IRDCs from Vasyunina et al. (2011) qualitatively agrees. Assuming uncertainties of about one order of magnitude for each of the derived abundances, only 2 of 9 molecules show a stronger deviation, namely, HNCO and SiO with lower abundances of slightly more than one order of magnitude in our sample. We detected CH3CN only in two sources in our IRDC sample, which was not detected by Vasyunina et al. either. In a more recent study, Vasyunina et al. (2014) additionally derived abundances of H2CO and CH3OH for a sample of IRDCs. The range of abundances they found are overlapping with, but slightly higher than, our findings. The observed sample of clumps in IRDCs from Sanhueza et al. (2012) shows in general higher abundances for all observed molecules, but 5 out of 6 of them agree with our values. Sakai et al. (2008) derived for a sample of IRDCs column densities of N2H+ with values that are similar to our results. For CH3OH our findings fall into the lower end of their observed range of column densities, but are still consistent. Sakai et al. (2010) calculated the column densities for HMPO-like sources. Comparing their sample of MSX sources with our HMPO sample we find good agreement with their observed range of column densities and our median values. Only CH3OH and CS isotopologues show a bigger deviation on the order of 10, which still agrees with the uncertainties of the observed values in our study. Our derived median abundances of the detections for the HMCs and UCHiis are similar to each other within the uncertainties. However, we derived higher median column densities for the HMCs due to a smaller mean distance and higher H2 column densities. Hatchell et al. (1998) inferred the H2 column densities from virial mass estimates of the observed gas, assuming a specific source size and extension. Therefore, their H2 column densities are different for different molecules. Comparing our median column densities for the HMCs and UCHiis with their results, we find an overall good agreement. The only molecule that deviates by a factor of more than 10 is CH3OH (for both phases) and CH3OCHO (for the UCHii stage). Reiter et al. (2011) derived molecular column densities of high-mass clumps coincident with water maser emission and mainly associated with Hii regions for a larger set of molecules. These values we can compare with our more evolved sources. The range of values they found for HCO+, C2H, HNC and CS are slightly lower, by a factor up to 5, than the median values in our sample. The column densities of N2H+ and SO are similar to our findings. The most prominent source in our sample is the well-studied Orion-KL hot core. A summary of abundances in this source is given in van Dishoeck & Blake (1998). Our results also agree well with these values. The SiO abundances differ by a factor of 6, while all other eight molecules in comparison show differences by factors of 14.

6. Conclusion

We observed the chemistry in 59 sources covering the early stages in massive star formation to study and characterize the physical and chemical evolution. This was the first consistent attempt to follow the evolution of physics and chemistry over several evolutionary steps in high-mass star formation. The underlying part of the evolutionary sequence assumed in this work started with an IRDC that evolved over the HMPO and HMC phase and ended with the formation of an UCHii region. Using this approach, we obtained a consistent sample of spectra for all of the sources. The spectra show an increasing number and intensity of the molecular transitions for more evolved stages. More complex and heavy molecules are formed with evolving age until the richness of molecular composition reaches its maximum in the HMC phase and declines for the UCHii stage, because these molecules are probably destroyed by the UV-radiation from the forming stars.

We calculated column densities for different molecular species that are in general increasing with evolutionary phase. Their ratios are less dependent on the temperature than the abundances. Thus we see a more obvious evolution in the abundances than in the ratios of these data. Furthermore, the derived lower limits of the ionization degree in the different stages increase along the evolutionary picture, as well.

In a second step we compared the characteristic median column densities of each evolutionary stage with a model. We calculated the best-fit parameters and ages, using a 1D physical model coupled to a gas-grain chemical model, in which the evolutionary sequence was approximately taken into account. The observed column densities and abundances agree reasonably well with the overall fit of our model. Furthermore, our results agree with previous observational studies. Based on the fit of the chemical model to 14 different molecular species of unresolved data, we were able to derive constrains on the underlying physical structure. The best-fit duration for the full evolution from an IRDC with atomic gas composition to the stage of an UCHii region is about 125 000 years. This is consistent with estimates from theoretical models of the formation of high-mass stars. The results from this study can be considered as chemical templates for the different evolutionary stages during high-mass star formation.


Acknowledgments

The authors wish to thank the anonymous referee for carefully reading the manuscript; his or her comments and suggestions improved it significantly. T.G. is supported by the Sonderforschungsbereich SFB 881 “The Milky Way System” (subproject B3) of the German Research Foundation (DFG). T.G. is member of the IMPRS for Astronomy & Cosmic Physics at the University of Heidelberg. This research made use of NASA’s Astrophysics Data System. D.S. acknowledges support by the Deutsche Forschungsgemeinschaft through SPP 1385: “The first ten million years of the solar system - a planetary materials approach” (SE 1962/1-2). T.V. wishes to thank the National Science Foundation (US) for its funding of the astrochemistry program at the University of Virginia. In the reduction and analysis of the data we made use of the GILDAS software public available at http://www.iram.fr/IRAMFR/GILDAS.

References

Online material

Table 1

Source list showing the position, the distance, and the evolutionary stage of all high-mass star-forming regions.

Appendix A: Appendix material

Table A.1

Integrated intensity in K km s-1.

Table A.2

Continuation of Table A.1 with additional species.

Table A.3

Iteration 0 column densities in a(x) = a × 10x cm-2 derived with the initially chosen typical temperatures (see Sect. 4.3 and Tables 36).

Table A.4

Continuation of Table A.3 with additional species.

Table A.5

Iteration 1 column densities in a(x) = a × 10x cm-2 derived with the mean temperatures from the best-fit models of iteration 0 (see Sect. 4.3 and Tables 36).

Table A.6

Continuation of Table A.5 with additional species.

thumbnail Fig. A.1

Modeled and observed median column densities in cm-2 in the IRDC stage. The modeled values are shown in black and the observed values in blue. The error bars are given by the vertical marks. The molecule name is labeled in the plot.

Open with DEXTER

thumbnail Fig. A.2

Modeled and observed median column densities in cm-2 in the HMPO stage. The modeled values are shown in black and the observed values in blue. The error bars are given by the vertical marks. The molecule name is labeled in the plot.

Open with DEXTER

thumbnail Fig. A.3

Modeled and observed median column densities in cm-2 in the HMC stage. The modeled values are shown in black and the observed values in blue. The error bars are given by the vertical marks. The molecule name is labeled in the plot.

Open with DEXTER

thumbnail Fig. A.4

Modeled and observed median column densities in cm-2 in the UCHii stage. The modeled values are shown in black and the observed values in blue. The error bars are given by the vertical marks. The molecule name is labeled in the plot.

Open with DEXTER

All Tables

Table 2

List of analyzed molecules with transitions, frequencies, energies of the upper level, critical densities, and effective densities calculated for 10 K and 100 Ka.

Table 3

Observed median abundances and standard deviation in a(x) = a × 10x for the IRDCs for T = 15 K and T = 20.9 K. x(e) is the ionization fraction.

Table 4

Observed median abundances and standard deviation in a(x) = a × 10x for the HMPOs for T = 50 K and T = 29.5 K. x(e) is the ionization fraction.

Table 5

Observed median abundances and standard deviation in a(x) = a × 10x for the HMCs for T = 100 K and T = 40.2 K. x(e) is the ionization fraction.

Table 6

Observed median abundances and standard deviation in a(x) = a × 10x for the UCHiis for T = 100 K and T = 36.0 K. x(e) is the ionization fraction.

Table 7

Initial atomic and molecular abundances.

Table 8

Parameters of the best-fit IRDC model.

Table 9

Parameters of the best-fit HMPO model.

Table 10

Parameters of the best-fit HMC model.

Table 11

Parameters of the best-fit UCHii model.

Table 12

Median column densities in a(x) = a × 10x for observations and best-fit IRDC model.

Table 13

Median column densities in a(x) = a × 10x for observations and best-fit HMPO model.

Table 14

Median column densities in a(x) = a × 10x for observations and best-fit HMC model.

Table 15

Median column densities in a(x) = a × 10x for observations and best-fit UCHii model.

Table 1

Source list showing the position, the distance, and the evolutionary stage of all high-mass star-forming regions.

Table A.1

Integrated intensity in K km s-1.

Table A.2

Continuation of Table A.1 with additional species.

Table A.3

Iteration 0 column densities in a(x) = a × 10x cm-2 derived with the initially chosen typical temperatures (see Sect. 4.3 and Tables 36).

Table A.4

Continuation of Table A.3 with additional species.

Table A.5

Iteration 1 column densities in a(x) = a × 10x cm-2 derived with the mean temperatures from the best-fit models of iteration 0 (see Sect. 4.3 and Tables 36).

Table A.6

Continuation of Table A.5 with additional species.

All Figures

thumbnail Fig. 1

From left to right, we show sample 16 GHz spectra of the sources IRDC048.6, HMPO18151, HMC029.96, and UCH013.87, representing the four selected evolutionary stages of high-mass star formation. The frequency resolution is ~0.2 MHz and typical rms values are ~0.1 K at 1 mm and ~0.03 K at 3 mm.

Open with DEXTER
In the text
thumbnail Fig. 2

Relative detection fractions for each analyzed molecular transition ordered from left to right: IRDCs without associated point sources at wavelengths below 70 μm, IRDCs with associated point sources, HMPOs, HMCs, and UCHiis. The color notation for the different source types is explained in the upper right corner.

Open with DEXTER
In the text
thumbnail Fig. 3

From top to bottom abundances (with respect to the total H2) of the analyzed molecules in the IRDC, HMPO, HMC and UCHii sample. The red line shows the median, the green cross is the mean, the bar indicates the inner 25%75% range around the median and the whiskers mark the total range of all detections. The black arrows indicate the lowest upper limit of all calculated upper limits for this particular molecule and stage.

Open with DEXTER
In the text
thumbnail Fig. 4

Snapshot of the radial distributions of CO, HCO+, N2H+, C2H, H2CO, CH3OH, CH3CN, and SiO at the time of the best-fit model of the HMC phase. The time is given relative to the beginning of the HMC phase.

Open with DEXTER
In the text
thumbnail Fig. 5

Relative abundances to H2 plotted for all four stages. The modeled values are shown as a black solid line, the observed values are depicted by a blue dashed line. The error bars are indicated by the vertical marks.

Open with DEXTER
In the text
thumbnail Fig. 6

Column density ratios for 15 different combinations of molecules plotted for all four stages. The modeled values are shown by the black solid line, the observed values are depicted by the blue dashed line. The error bars are indicated by the vertical marks.

Open with DEXTER
In the text
thumbnail Fig. 7

Averaged relative abundances to H2 for electrons, HCO+, and N2H+ plotted for the all four stages. The observed values are shown as solid lines, the modeled values are depicted as dashed lines. The electron abundances are plotted in red, N2H+ abundances in black, and HCO+ abundances in blue. The error bars for both the model and the observational data are a factor of about 3 (not shown).

Open with DEXTER
In the text
thumbnail Fig. A.1

Modeled and observed median column densities in cm-2 in the IRDC stage. The modeled values are shown in black and the observed values in blue. The error bars are given by the vertical marks. The molecule name is labeled in the plot.

Open with DEXTER
In the text
thumbnail Fig. A.2

Modeled and observed median column densities in cm-2 in the HMPO stage. The modeled values are shown in black and the observed values in blue. The error bars are given by the vertical marks. The molecule name is labeled in the plot.

Open with DEXTER
In the text
thumbnail Fig. A.3

Modeled and observed median column densities in cm-2 in the HMC stage. The modeled values are shown in black and the observed values in blue. The error bars are given by the vertical marks. The molecule name is labeled in the plot.

Open with DEXTER
In the text
thumbnail Fig. A.4

Modeled and observed median column densities in cm-2 in the UCHii stage. The modeled values are shown in black and the observed values in blue. The error bars are given by the vertical marks. The molecule name is labeled in the plot.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.