© ESO, 2013

1. Introduction

Millimeter lines from complex organic molecules are widely associated with high-mass star forming regions and indeed form one of the signposts of the deeply embedded phase of star formation (e.g., Blake et al. 1987; Hatchell et al. 1998; Gibb et al. 2000; Fontani et al. 2007; Requena-Torres et al. 2008; Belloche et al. 2009; Zernickel et al. 2012). Many studies of the chemistry in such regions have been carried out, either through complete spectral surveys of individual sources or by targetting individual molecules in a larger number of sources (see Herbst & van Dishoeck 2009; Caselli & Ceccarelli 2012, for reviews). In spite of all this work, only few systematic studies of the abundances of commonly observed complex molecules have been performed across a sample of massive young stellar objects (YSOs), to search for similarities or differences depending on physical structure and evolutionary state of the object. Intercomparison of published data sets is often complicated by the use of different telescopes with different beams, different frequency ranges and different analysis techniques.

Chemical abundances depend on the physical structure of the source such as temperature, density and their evolution with time, as well as the amount of UV radiation impinging on the gas and dust. In contrast with the case of solar-mass stars, the physical structures and mechanisms for forming massive (M > 8 M_⊙) stars are still poorly understood. Indeed, theoretically, the powerful UV radiation pressure from a high-mass protostellar object (HMPO) should prevent further accretion and so inhibit the formation of more massive stars (Zinnecker & Yorke 2007). However, a number of recent studies have claimed the presence of disk- or toroid-like “equatorial” structures surrounding a handful of high-mass protostars (Cesaroni et al. 2007). These data support theories in which high-mass star formation proceeds in a similar way as that of low-mass stars via a disk accretion phase in which high accretion rates and non spherically symmetric structures overcome the problem of radiation pressure. The best evidence so far is that for a ~5000 AU disk in Keplerian rotation around IRAS 20126+4104, claimed by Cesaroni et al. (1999) on the basis of the presence of a velocity gradient in the CH₃CN emission perpendicular to the direction of the outflow, as predicted by the disk-accretion paradigm. Surprisingly, even in the best case of IRAS 20126+4104, the detailed chemistry of these (potential) disks has not yet been studied.

Most chemical models invoke grain surface chemistry to create different generations of complex organic molecules (Tielens & Charnley 1997). Hydrogenation of solid O, C, N and CO during the cold (T_d < 20 K) prestellar phase leads to ample production of CH₃OH and other hydrogenated species like H₂O, CH₄ and NH₃ (Tielens & Hagen 1982). Exposure to UV radiation results in photodissociation of these simple ices, with the fragments becoming mobile as the cloud core heats up during the protostellar phase. First generation complex molecules result from the subsequent recombination of the photofragments, and will eventually evaporate once the grain temperature rises above the ice sublimation temperature of ~100 K (Garrod & Herbst 2006; Garrod et al. 2008). Good examples are C₂H₅OH, HCOOCH₃ and CH₃OCH₃ resulting from mild UV processing of CH₃OH ice (Öberg et al. 2009). Finally, a hot core gas phase chemistry between evaporated molecules can drive further complexity in second generation species (e.g., Millar et al. 1991; Charnley et al. 1992, 1995).

One of the most obvious consequences of an equatorial rather than spherical structure is that UV radiation can more easily escape the central source and illuminate the surface layers of the surrounding disk or toroid as well as the larger scale envelope (Bruderer et al. 2009, 2010) (Fig. 1). This can trigger enhanced formation of complex organic molecules in the ices relative to methanol. Another effect of UV radiation is that it drives increased photodissociation of gaseous N₂ and CO. The resulting atomic N and C would then be available for grain surface chemistry potentially leading to enhanced abundances of species like HNCO and NH₂CHO.

To investigate the effects of disk-like structures on the chemistry, we present here a single-dish survey using the James Clerk Maxwell Telescope (JCMT) of three HMPOs for which large equatorial structures (size > 2000 AU) have been inferred, namely IRAS 20126+4104, IRAS 18089–1732, and G31.41+0.31. The results are compared with those of a recent survey of a sample of HMPOs by Bisschop et al. (2007, hereafter BIS07), targeting many of the same molecules and settings. The use of the same telescope and analysis method allows a meaningful comparison between the two samples of sources. BIS07 found that the O-rich complex molecules are closely correlated with the grain surface product CH₃OH supporting the above general chemical scenario. N-rich organic molecules do not appear to be correlated with O-rich ones, but overall, the relative abundances of the various species are found to be remarkably constant within one chemical family. One of the main questions is whether this similarity in abundance ratios also holds for sources with disk-like structures. Although our observations do not spatially resolve these structures, they are sensitive enough and span a large enough energy range to determine abundances on scales of ~1″ and thereby set the scene for future high-angular resolution observations with interferometers like the Atacama Large Millimeter/submillimeter Array (ALMA). Moreover, current interferometers resolve out part of the emission, which is why single-dish observations are still meaningful.

This paper is organized as follows. In , the observed sources and frequency settings are introduced and the details of the observations are presented. focuses on the data analysis methods. Specifically, two techniques are used to determine excitation temperatures and column densities: the rotation diagram method employed by BIS07 and spectral modeling tools in which the observed spectra are simulated directly. presents the results from both analysis methods and discusses their advantages and disadvantages. Finally, compares our results with those of BIS07 to draw conclusions on similarities and differences in chemical abundances between sources with and without large equatorial structures.

2. Observations

2.1. Observed sources

Table 1 gives the coordinates, luminosity L, distance d, galactocentric radius R_GC, velocity of the local standard of rest V_LSR and the typical line width for the observed sources. The selected sources are massive YSOs, for which strong evidence exists for a circumstellar disk structure. All sources are expected to harbor a hot core: a compact, dense (≥10⁷ cm^-3) and warm (≥100 K) region with complex chemistry triggered by the grain mantle evaporation (Kurtz et al. 2000). CH₃CN emission from ≥100 K gas is present in all sources. Moreover, CH₃OH 7_K-6_K transitions (338.5 GHz) with main beam temperatures of ≥1 K have been observed for these sources. Sources also needed to be visible from the JCMT¹.

2.1.1. IRAS 20126+4104

IRAS 20126+4104 (hereafter IRAS 20126) is a luminous (~10⁴L_⊙) YSO located relatively nearby at a distance of 1.64 kpc (Moscadelli et al. 2011). It was first identified in the IRAS point source catalog by the IR colors typical of ultracompact HII regions and by H₂O maser emission characteristic of high-mass star formation (Comoretto et al. 1990). IRAS 20126 features a ~0.25-pc scale inner jet traced by H₂O maser spots in the SE-NW direction with decreasing velocity gradient (Tofani et al. 1995). Source and masers are embedded inside a dense, parsec-scale molecular clump (Estalella et al. 1993; Cesaroni et al. 1999). The inner jet feeds into a larger scale bipolar outflow with the two having reversed velocity lobes (Wilking et al. 1990; Cesaroni et al. 1999). The reversal is likely to be due to precession of the inner jet caused by a companion separated by a distance of ~0.5″ (850 AU) (Hofner et al. 1999; Shepherd et al. 2000; Cesaroni et al. 2005; Sridharan et al. 2005). A rotating, flattened, Keplerian disk structure has been detected perpendicular to the inner jet. Observations of CH₃CN transitions show a Keplerian circumstellar disk (radius ~1000 AU) with a velocity gradient perpendicular to the jet and a hot core with a diameter of ~0.0082 pc and a temperature of ~200 K at a geometric center of the outflow (Cesaroni et al. 1997, 1999; Zhang et al. 1998). Direct near-infrared (NIR) observations show a disk structure as a dark line (Sridharan et al. 2005). The disk shows a temperature and density gradient and is going through infall of material with a rate of ~2 × 10^-3M_⊙  yr^-1 as expected for a protostar of this mass and luminosity (Cesaroni et al. 2005). Recent modeling of the Spectral Energy Distribution (SED) of IRAS 20126 shows indeed a better fit when a disk is included (Johnston et al. 2011).

2.1.2. IRAS 18089-1732

IRAS 18089-1732 (hereafter IRAS 18089) is a luminous (~10^4.5L_⊙, Sridharan et al. 2002) YSO located at a distance of 2.34 kpc (Xu et al. 2011). It was identified based on CS detections of bright IRAS point sources with colors similar to ultracompact HII regions and the absence of significant free-free emission (Sridharan et al. 2002), and with H₂O and CH₃OH maser emission (Beuther et al. 2002). The CO line profile shows a wing structure characteristic of an outflow, although no clear outflow structure could be resolved from the CO maps (Beuther et al. 2002). A collimated outflow in the Northern direction is however seen in SiO emission on scales of 5″ (Beuther et al. 2004). IRAS 18089 has a highly non-circular dust core of ~3000 AU diameter (~1″), with optically thick CH₃CN at ~350 K (Beuther et al. 2005). HCOOCH₃ was found to be optically thin, with emission confined to the core, and showing a velocity gradient perpendicular to the outflow indicative of a rotating disk (Beuther et al. 2004). Also hot NH₃ shows a velocity gradient perpendicular to the outflow, although no Keplerian rotation was found, possibly due to infall and/or self gravitation (Beuther & Walsh 2008). Several hot-core molecules (HCOOCH₃, CH₃CN, CH₃OCH₃, HNCO, NH₂CHO, CH₃OH, C₂H₅OH) were mapped by Beuther et al. (2005) but no column densities or abundances were reported.

2.1.3. G31.41+0.31

G31.41+0.31 (hereafter G31) is a luminous (2.6 × 10⁵L_⊙, Cesaroni et al. 1994b) YSO at a distance of 7.9 kpc (Churchwell et al. 1990). Preliminary evidence for a rotating disk with a perpendicular bipolar outflow was reported by Cesaroni et al. (1994a,b) and Olmi et al. (1996) showing a velocity gradient across the core in the NE-SW direction, similar to previously detected OH masers (Gaume & Mutel 1987). High-angular resolution CH₃CN observations by Beltrán et al. (2005) could not detect Keplerian rotation typical for less luminous stars, nevertheless a toroidal structure undergoing gravitational collapse and fast accretion (~3 × 10^-2M_⊙ yr^-1) onto the central object was found. The G31 hot molecular core (HMC) is part of a complex region in which multiple stellar sources are detected (Benjamin et al. 2003); indeed, it is separated from an ultracompact (UC) HII region by only ~5″ and overlaps with a diffuse halo of free-free emission, possibly associated with the UC HII region itself (Cesaroni et al. 1998). More recent interferometric observations confirm the velocity gradient in the NH₃ (4,4) inversion transition and in CH₃CN data (Cesaroni et al. 2010, 2011). Line profiles look like a rotating toroid with infall motion. Several complex hot-core molecules have been observed in G31, including glycolaldehyde CH₂OHCHO (Beltrán et al. 2009), but again no abundance ratios have been presented.

2.2. Observational details

The observations were performed at the JCMT on Mauna Kea, Hawaii, between August 2007 and September 2009. The observations of the 338 GHz region covering CH₃OH (7_K → 6_K) transitions were taken from JCMT archive. The front ends consisted of the facility receivers A3 (230 GHz region) and HARP-B (340 GHz region). The back-end was the digital autocorrelation spectrometer (ACSIS), covering 400 and 250 MHz of instantaneous bandwidth for A3 and HARP-B, respectively, with a channel width of 50 kHz. The noise level for both receivers was T_rms ~ 20 mK on a scale when binned to 0.5 km s^-1. The integration time was ~1 h and 1.8 h for A3 and HARP-B, respectively. The spectra were scaled from the observed antenna temperature, , to main-beam temperature, T_MB, using main beam efficiencies η_MB of 0.69 and 0.63 at 230 GHz and 345 GHz, respectively. We adopt a calibration error of 20%.

The HPBW (half-power beam width, θ_B) for the 230 and 345 GHz band observations are 20–21″ and 14″, respectively. Emission from a volume with a source diameter θ_S ≤ θ_B undergoes beam dilution described by the beam-filling factor, η_BF: (1)Table 2 gives the observed frequency settings and the targeted molecular lines. The settings were taken from BIS07 and were chosen to cover at least one strong line for the target molecules as well as lines of other interesting species. Strong lines of target molecules were chosen due to their high main-beam temperatures and minimum line-confusion in line surveys of Orion-KL by Sutton et al. (1985) and Schilke et al. (1997) at 230 GHz and 345 GHz, respectively. In order to determine rotation temperatures, we covered at least two transitions for a given species with E_up < 100 K and E_up > 100 K each. BIS07 used the single pixel receiver B3, the predecessor of HARP-B, together with DAS (Digital Autocorrelation Spectrometer) as the back-end, covering a larger instantaneous bandwidth of 500 MHz. Only the central receptor of HARP-B array is analyzed here as no significant off source emission was detected in the complex molecules.

Data reduction and line fitting were done using the CLASS software package². Line assignments were done by comparison of observed frequencies corrected for source velocity with the JPL³, CDMS⁴ and NIST⁵ catalogs (Pickett et al. 1998; Müller et al. 2005). The line assignment/detection was based on Gaussian fitting with the following criteria: (i) the fitted line position had to be within ±1 MHz of the catalog frequency; (ii) the FWHM had to be consistent with those given in Table 1; and (iii) a S/N > 3 is required on the peak intensity. Section A.1 in the Appendix describes in more detail the error analysis. In general, our errors on the integrated intensities are conservative and suggest a lower S/N than that on the peak intensity or obtained using more traditional error estimates.

3. Data analysis

3.1. Rotation diagrams

Rotation temperatures and column densities were obtained through the rotation diagram (RTD) method (Goldsmith & Langer 1999), when 3 or more lines are detected over a sufficiently large energy range. Integrated main-beam temperatures, , can be related to the column density in the upper energy level by: (2)where N_up is the column density in the upper level, g_up is the degeneracy of the upper level, k is Boltzmann’s constant, ν is the transition frequency, μ is the dipole moment and S is the line strength. The total beam-averaged column density N_T in cm^-2 can then be computed from: (3)where Q(T_rot) is the rotational partition function, and E_up is the upper level energy in K.

Blended transitions of a given species with similar E_up (ΔE_up < 30 K) were assigned intensities according to their Einstein coefficients (A) and upper level degeneracies (g_up): (4)where the summation is over all the contributing transitions. Blended transitions with different E_up, or with contamination from transitions belonging to other species, were excluded from the RTD fit.

The beam averaged column density, N_T, is converted to the source-averaged column density N_S using the beam-filling factor η_BF: (5)The standard RTD method assumes that the lines are optically thin. Lines with strong optical depth, determined from the arguments in Sect. 4.1 as well as the models discussed in Sect. 3.2, were excluded from the fit. For CH₃OH, also low-E_up lines arising from a cold extended component (see Sect. 4.4.1) were excluded.

Differential beam dilution is taken into account by multiplying the line intensities in the 230 GHz range by η_BF (340 GHz)/ η_BF(230 GHz). For warm and cold molecules, beam dilution is derived assuming source diameters θ_{T = 100  K} (see Eq. (10) below) and 14″, respectively (see Eq. (1)). All emission is thus assumed to be contained within the smallest beam size. Same approach was used in BIS07.

The vibrational partition function was ignored assuming that all emission arises from the ground vibrational state so that the vibrational partition function can be set to unity. This approximation can cause an underestimation of the derived column densities for larger molecules, even at temperatures of 100–200 K. Indeed, Widicus Weaver et al. (2005) show that the error can be up to a factor of 2 for temperatures up to 300 K. Since this approximation affects all complex molecules (albeit at different levels), the overall error in abundance ratios will be less than a factor of 2 and well within the other uncertainties.

3.2. Spectral modeling

The alternative method for analysing the emission is to model the observed spectra directly. For this purpose, we used the so-called “Weeds” extension of the CLASS software package⁶, developed to facilitate the analysis of large millimeter and submillimeter spectral surveys (Maret et al. 2011). In this model, the excitation of a species is assumed to be in LTE (Local Thermodynamic Equilibrium) at a temperature T_ex. The brightness temperature, T_B, of a given species as a function of the rest frequency ν is then given by: (6)where η_BF is the beam filling factor (see Eq. (1)), J_ν is the radiation field such that:

and T_bg is the temperature of the background emission. The HPBW θ_B is calculated within the Weeds model as θ_B = 1.22c/νD, where c is the speed of light and D is the diameter of the telescope⁷. The opacity τ(ν) is: (7)where the summation is over each line i of the considered species. is the rest frequency of the line and φⁱ is the profile function of the line, given by: (8)where Δν is the line width in frequency units at 1/e height. Δν can be expressed as a function of the line FWHM in velocity units ΔV by: (9)The input parameters in the model are the column density N_S in cm^-2, excitation temperature T_ex in K, the source diameter θ_S in arcseconds, offset velocity from the source LSR (Local Standard of Rest) in km s^-1 and the line FWHM in km s^-1. All parameters excluding the source diameter are optimized manually to obtain the best agreement with the observed spectra. The source diameter for emission from cold species is allowed to vary over a larger area; generally 14″ is used. The emission from hot-core molecules targeted in our work is assumed to originate from the central region with T_dust ≥ 100 K. The source diameter for the warm emission is calculated using a relation derived from dust modeling of a large range of sources (BIS07): (10)where L/L_⊙ is the luminosity of the source relative to the solar luminosity. Table 3 gives the calculated R_{T = 100  K} radii and diameters for the observed sources.

In the analysis for individual molecules, the initial values for T_ex were based on the T_rot from the RTD analysis in case of optically thin species. For optically thick species the value of T_ex from the ¹³C-isotope was used. If T_ex for the isotopologue could not be obtained, the initial temperature was guessed. The T_ex value was then optimized visually based on the relative intensity of the emission lines. The simulated emission was not allowed to exceed the emission of optically thin lines in the spectrum in any of the observed spectral ranges. Coinciding and blended transitions, which together contribute to an optically thick line, are excluded in the analysis.

Due to the visual optimization and the possibility of overlapping lines (particularly in the line-rich source G31), the resulting T_ex values are only a rough estimate (±50–100 K) and do not differ significantly from those from the RTD method. The column density for each species was constrained by optically thin, unblended lines, where available.

For the specific case of CH₃OH, which requires two temperature components for a proper fit, the analysis was also done with the CASSIS analysis package⁸. CASSIS has the advantage that it can properly model the emission from overlapping optically thick lines, as well as from nested regions of emission. For CH₃OH, warm emission from the compact inner region may be absorbed by the surrounding colder gas, which can influence the derived hot core abundances.

4. Results

4.1. General results and comparison between sources

The observed sources, IRAS 20126, IRAS 18089 and G31, differ from each other significantly in the observed chemical richness in the JCMT single-dish data. Figure 2 presents two frequency ranges with lines from several observed species. For G31, strong lines of all complex organic species are detected, whereas for IRAS 20126 many targeted lines are below the detection limit. Many complex molecules are also found in IRAS 18089 but with weaker lines than for G31. Among the serendipitous discoveries acetone, CH₃COCH₃, is possibly seen in G31 and IRAS 18089 for the first time (see Sect. 4.5). Integrated line intensities for all detected lines and selected upper limits are given in Tables A.1–A.14. The rotation diagrams for the detected species are shown in Figs. B.1–B.13 whereas the optimized parameters in the Weeds model for each molecule and source are available online in Tables C.1–C.3.

4.2. Optical depth determinations

To assess the importance of optical depths effects, the ratio of lines of different isotopologues are compared. The expected ¹²C/¹³C isotope ratio depends on the galactocentric radius, R_GC, according to Eq. (11) (Wilson & Rood 1994). (11)The isotope ratios derived for our sources using Eq. (11) are given in Table 1. The galactocentric radii were calculated trigonometrically from the galactic coordinates, using the IAU value for the distance to the Galactic center R₀ = 8.5 kpc (Kerr & Lynden-Bell 1986).

Table 4 lists the observed isotopologue intensity ratios for the most abundant species in our sources. The CH₃OH/¹³CH₃OH ratios are derived from a low-energy transition 2_2,0,+0–3_1,3,+0 with E_up = 44.6 K. No high E_up transitions were reliably detected for ¹³CH₃OH in our sources in our standard settings and the low E_up ratios are therefore taken to apply to cold methanol. For G31.41+0.31, one additional observation was carried out to cover a transition with E_up > 100 K. A CH₃OH/¹³CH₃OH intensity ratio of 4.0 is derived for a transition with E_up ≈ 210 K (13_1,12,−0–13_0,13,+0). This indicates that also warm methanol is optically thick in this source.

The CH₃CN/CH₃¹³CN ratios are derived from the 13₃–12₃ line intensities for G31 and IRAS 18089 and indicate that CH₃CN is optically thick in these two sources, but not in IRAS 20126. For H₂CO and HNCO, isotopologue lines are detected but for different transitions than the main isotopologues. Thus, a model is needed to infer the optical depths. Fits to each of the isotopologues independently at a fixed temperature of 150 K using the R_{T = 100  K} source size gives column density ratios that are significantly smaller than the overall isotope ratios, suggesting that these species are also optically thick for G31 and IRAS 18089. In practice, we have excluded the optically thick lines (as indicated by the Weeds model) from the RTD fits for all species.

4.3. Temperatures

Table 5 summarizes the derived temperatures from the RTD fit for the various species. As also found by BIS07, molecules can be classified into cold (<100 K) and warm (>100 K), and our categorization is similar to theirs. The Weeds analysis is consistent with this grouping. There is however variation in temperatures within the groups, with warm species having rotation temperatures from 70 to 300 K, and cold species from 40 to 100 K. Some variation is seen in rotation temperatures of individual species between different sources; the rotation temperatures are generally higher in G31 than in IRAS 18089, while IRAS 20126 has the lowest of the three.

The T_rot value for CH₃OH is ~300 K for G31 and IRAS 18089 and ~100 K for IRAS 20126. For G31, several lines with E_up > 400 K are detected, which makes the RTD fit more robust. For IRAS 18089 and IRAS 20126 the accuracy of the RTD fits suffers from the small range of E_up in the detected transitions. In addition to optically thick lines, low-E_up lines have been excluded from the fit. These lines are underestimated by the RTD fit and probably originate from a colder extended region also seen in the ¹³C lines (Fig. 3). The T_rot from the RTD analysis therefore represents the warm CH₃OH alone. See Sect. 4.4.1 for a more detailed discussion on the CH₃OH emission.

For CH₃CN, the T_rot values range from ~200 and ~350 K for the three sources. The CH₃¹³CN rotation diagram gives a value of T_rot of only 53 K, however (Fig. 3). This illustrates the large uncertainties at high temperatures for optically thick species and the possibility of a cold component in addition to the warm one.

Contrary to the general trend, the T_rot value for H₂CO is somewhat higher (204 K) in IRAS 18089 than in G31 (157 K). The discrepancy could be influenced by the small number of lines used for the analysis. All fitted lines are those belonging to the para-H₂CO, and the T_rot fits are thus not affected by fluxes from different spin states. A single transition of H₂¹³CO (3_1,2–2_1,1) was covered, and no information on the excitation temperature can therefore be obtained from the minor isotopologue.

The T_rot values of the other complex species, HNCO, C₂H₅OH, C₂H₅CN, NH₂CHO and CH₃OCH₃ are around 100 K and are slightly higher in G31 than in IRAS 18089. The T_rot of HCOOCH₃ stands out in both sources, in G31 being closer to 200 K. No lines belonging to these species were detected in IRAS 20126.

The species classified as cold by BIS07, CH₂CO and CH₃CCH, indeed show cold rotation temperatures in all sources. Not enough lines of CH₃CHO or HCOOH, which are also classified as cold in BIS07, are observed in our sources for making rotation diagrams.

4.4. Column densities

Table 7 presents the column densities derived from the RTD analysis, Weeds or CASSIS model, and from the ¹³C-isotopologues for the optically thick species. Following BIS07, the column densities for warm molecules are given as source-averaged values (see Eq. (5)). The emission from cold molecules extends over a larger volume and the values are given as beam-averaged column densities. Typical uncertainty of the column densities derived from the RTD analysis is ~40%.

4.4.1. CH₃OH

An accurate determination of the CH₃OH column density is essential for comparing the abundance ratios of complex organic species. For hot-core molecules, it is particularly important to quantify the warm CH₃OH emission. The column densities of CH₃OH in BIS07 were determined by the RTD method excluding the optically thick lines. The same is done in our analysis. Our rotation diagrams however also show emission from low-E_up transitions, which are strongly underestimated by the RTD fit on the warm lines, providing further evidence for the presence of a colder component. We have therefore also excluded these transitions. The fit to the higher E_up, optically thin lines should give the warm CH₃OH column density obtained in a similar way as BIS07.

CH₃OH emission was also simulated using a two-component CASSIS model. A single-component model is not able to simultaneously reproduce the warm and cold lines without overestimating the lines from intermediate energy levels. Indeed, a better agreement is obtained using a two-component model, consisting of a warm compact component and a cold extended component. Table 6 shows the best model parameters obtained from the fitting. The warm components are fixed to θ_{T = 100  K} while the cold component is allowed to extend beyond the beam diameter. The warm CH₃OH column densities derived from the CASSIS fit are in agreement with the values derived from the RTD analysis. The best fits plotted onto the CH₃OH 7_K − 6_K transitions (338.5 GHz) are shown in Fig. 4.

Several ¹³CH₃OH lines are detected in G31. However, only low-E_up lines are reliably detected since the high-E_up lines are very weak or blended. Assigning these lines to the cold component and assuming T_rot = 20 K (Öberg et al. 2011; Requena-Torres et al. 2008), the beam-averaged ¹³CH₃OH column density is 2.1 × 10¹⁵ cm^-2, corresponding to a CH₃OH column density of 1.7 × 10¹⁷ cm^-2 for the 14″ volume. This is more than an order of magnitude higher than found in the Weeds model for the cold component, supporting the optically thick interpretation of the cold component. Similarly, for the case of IRAS 18089, the detected ¹³CH₃OH lines may indicate a high optical depth.

As mentioned before, the ¹³CH₃OH line with E_up = 211 K, covered in additional observations for G31, reveals optical thickness in the warm component as well. The CH₃OH column density derived from this line assuming temperature of 300 K is 2.7 ×   10¹⁸ cm^-2 and the CASSIS fit as well as the RTD analysis underestimate the warm column density by a factor of ~2.5. In principle, the CASSIS analysis could be made consistent with the ¹³CH₃OH value by letting the warm source size vary as well. However, since such information is not known for other molecules, as well as for consistency with BIS07, we have chosen to keep the warm source size fixed at the 100 K radius.

In summary, it seems plausible that at least in our sources, methanol emission arises from two temperature components, a warm (T_ex ≈ 300 K) compact component and a cold (T_ex ≈ 20 K) significantly more extended component. At least for G31, the CH₃OH emission is optically thick in both components.

4.4.2. Other molecules

Overall, the column densities of the various species follow the same pattern in all sources, and hence seem to be well correlated relative to each other. CH₃OH is the most abundant complex molecule in all sources. The other species have in general half to one order of magnitude lower column densities.

The CH₃CN emission is optically thick for G31 and IRAS 18089. Column densities from the RTD analysis are therefore underestimated. Due to a lack of sufficient optically thin lines, also the Weeds analysis underestimates the column densities. Values derived from the ¹³C-isotopologues are thus more reliable, even though some fraction may arise from a colder component. Indeed, the column density for CH₃CN derived from the 13-isotopologue is an order of magnitude larger than that from the RTD of the main isotopologue alone. The same procedure was used by BIS07.

For H₂CO, the best estimates come from the RTD analysis on the optically thin para-H₂CO lines. The derived H₂CO column densities assuming T_rot from Table 5 and a statistical ortho-to-para ratio of 3 are given in Table 7. The column densities derived from the ¹³C-isotopologue of ortho-H₂CO are still larger than from RTD analysis using optically thin lines, which indicates either larger ortho-to-para ratio in these sources (ratio of ~3 to 5 has been predicted for cold clouds by Kahane et al. 1984) or non LTE excitation. BIS07 derived the H₂CO column densities from the ¹³C-isotopologue.

The column density derived for HCOOCH₃ from the RTD diagram is significantly smaller than that derived from the Weeds model. The HCOOCH₃ emission is stronger in the 230 GHz beam (20–21″) than in the 345 GHz beam (14″) probably indicating significant extended emission. The RTD analysis was performed on the entire dataset, while only lines in the 230 GHz band was used in the Weeds modeling.

5. Discussion

5.1. Comparison to massive YSOs without a disk structure

In order to see the effect of a flattened circumstellar structure on the chemistry around the YSO, we compare the temperatures and abundances of the complex organic molecules in sources studied in this work and those studied in BIS07. In addition to our sources (IRAS 20126, IRAS 18089 and G31), AFGL 2591, NGC 7538 IRS1 and G24.78 from the BIS07 sample are now known to have disk structures as well (van der Tak et al. 2006; Wang et al. 2012; Pestalozzi et al. 2004, 2009; Knez et al. 2009; Beltrán et al. 2005). A disk-like structure may be present also for other sources but current evidence is not as strong as for the above sources. For the comparison with the BIS07 sample, we use the results primarily from the RTD method to exclude method based differences.

5.1.1. Rotation temperatures

Figure 7 shows the rotation temperatures of the complex species for the observed sources together with those from BIS07. As in BIS07, we find that the complex molecules can be divided into warm and cold species based on their rotation temperatures. The rotation temperatures obtained for the molecules in our sources agree with the division; CH₃OH, H₂CO, C₂H₅OH, HCOOCH₃, CH₃OCH₃, HNCO, NH₂CHO, CH₃CN and C₂H₅CN show rotation temperatures generally ≥100 K while CH₂CO, CH₃CHO, HCOOH and CH₃CCH have rotation temperatures <100 K.

The rotation temperatures of the cold molecules show a very small scatter (±25 K) from source to source, with or without a disk structure, ignoring the one outlier for HCOOH. Of the warm species, C₂H₅OH has a consistent rotation temperature of 100–150 K from source to source. The RTD analysis for these species is reliable due to low optical depths and lack of anomalous excitation.

The other warm species show more scatter in the derived rotation temperatures. This is particularly the case for CH₃OH and CH₃CN, which have rotation temperatures ranging from 100 to 350 K. No systematic difference is seen between the two source types. Moreover, the HCOOCH₃ rotation temperatures for AFGL 2591 and NGC 7538 IRS1 are similar to the diskless sources, disproving any difference between the two source types. The scatter in the rotation temperatures could indicate that they exist in environments with different temperatures, but it may also be caused by optical depth effects (H₂CO, HNCO, CH₃OH and CH₃CN) and (anomalous) radiative excitation. The latter has been previously seen for, e.g., NH₂CHO and HCOOCH₃ (Nummelin et al. 2000) as well as HNCO (Churchwell et al. 1986). For species with significant radiative excitation, e.g., HCOOCH₃, the scatter of the data points in the RTD plots is higher and results in additional uncertainty since the linear fits do not properly capture this scatter.

5.1.2. Column densities

Figures 8 and 9 show the column densities for the targeted species in our sources and those from BIS07 for warm and cold molecules, respectively. For the warm species the column densities within the sources vary by 1–2 orders of magnitude. There is also a significant variation between sources. G31 is chemically richest of the sources, with highest column densities of all the targeted species, compared to other sources. IRAS 20126 is among the chemically poorest sources. The pattern of column densities is remarkably similar, however, and the disk sources do not stand out.

5.1.3. Abundance ratios

Table 8 gives the abundance ratios for the targeted species with respect to CH₃OH for oxygen-bearing species and HNCO for nitrogen-bearing species. Column densities from the RTD analysis are primarily used. In cases where column densities are not available from the RTD analysis (mainly for IRAS 20126) the values are taken from the Weeds analysis.

The resulting abundance ratios are presented in Fig. 10. Values from the RTD method are shown in black bars and can thus be directly compared to BIS07, whereas those from the Weeds model and ¹³C-isotopologue are in red and blue bars. In general, the abundance ratios in different sources have larger variations within the source types than between them. For example, H₂CO has both lowest and highest abundance ratio for two disk candidates, G31 and IRAS 20126, respectively. The H₂CO results depend somewhat on the analysis method used.

C₂H₅OH, HCOOCH₃ and CH₃OCH₃ abundances peak for G31 and IRAS 18089, and are generally lower for other sources. CH₃OCH₃ also peaks for some diskless sources. The N-bearing species NH₂CHO has a large variation in the abundance ratio with respect to HNCO. Among our sources, NH₂CHO peaks for G31, the source with largest abundance of O-bearing species, and the least clear disk structure.

Figure 11 shows the abundance ratios of N-bearing species with respect to CH₃OH. Again, a lot of variation is seen within each source type, with no specific trend found for sources with a disk-like structure.

We stress that the absolute inferred column densities and abundance ratios are uncertain by a factor of a few up to an order of magnitude, as already indicated by the different analysis techniques. An independent assessment of the accuracy of the results can be made by comparison with the inferred column densities and abundances of Zernickel et al. (2012) for NGC 6334 I (=NGC 6334 IRS), who observed a completely different set of lines of the same molecules with Herschel-HIFI and the SMA in beams ranging from 2–40″. Comparison with the results of BIS07 shows good agreement within a factor of a few for the abundances of several species (CH₃OH, CH₃OCH₃, CH₃OCHO, C₂H₅OH, C₂H₅CN) whereas others differ by an order of magnitude (H₂CO, CH₃CN). The species that show the largest discrepancy are those with highly optically thick lines and without a large set of isotopologue lines. Thus, a combination of differences in adopted source sizes and optical depth effects can account for the discrepancies. Because our approach is the same for all sources, the relative values from source to source should be more reliable.

5.2. Chemical and physical implications

The main goal of this study is to investigate similarities and differences between sources with and without a disk-like structures. The presence of a flattened accretion disk should allow the UV radiation from the central object to escape more readily and then impinge on the gas and dust in the outer walls of a flared disk or outflow cavity (Bruderer et al. 2010). Increased UV radiation could manifest itself as enhancement in the complex organic species produced through UV photodissociation from CH₃OH in the solid state, such as CH₃OCH₃, C₂H₅OH and CH₃OCHO (Öberg et al. 2009). In this scenario, CH₃OH is photodissociated into various radicals such as CH₃, CH₃O and CH₂OH which can become mobile at higher ice temperatures (20–40 K) and form the observed complex organic molecules (Garrod & Herbst 2006). Higher temperatures favor diffusion of larger radicals resulting in the formation of larger complex organic molecules compared with small molecules like H₂CO. Another related parameter that plays a role is the CO content of ices, with H₂CO, CHO- and COOH-containing species enriched in cold, CO-rich CH₃OH ices in which CO has not evaporated. Indeed Öberg et al. (2011) show that the increased abundance of HCOOCH₃ in low-mass YSOs compared with the high-mass sample of BIS07 may be due to this effect. Although only limited laboratory experiments exist on N-containing molecules, the abundance of N-bearing complex organic molecules could be enhanced due to photodissociation of N₂ or NH₃.

Our results show that there is no consistent enhancement of complex organic species produced through CH₃OH photochemistry, nor of N-bearing species, that can be attributed to the presence of a flattened disk structure on the scales probed by our data. The variations within source types are larger than between the source types. One explanation for this lack of differentation could be that all high-mass sources have a source structure through which UV radiation can impact the larger surroundings. In particular, just like low-mass sources, all high-mass YSOs are expected to have outflows and cavities through which UV radiation can escape and affect the chemistry, whether or not they have a large disk-like structure. This scenario can be tested with future high angular resolution data on <1″ scale with ALMA which should then reveal the emission from complex molecules coating the walls of the outflow cavities.

An alternative explanation is that enhanced temperature in the YSO environment is not needed for the production of complex organic molecules and that they are already formed in the prestellar stage through UV- and cosmic-ray processing of cold ices, using just the internal UV field produced by interaction of cosmic rays with H₂ (Prasad & Tarafdar 1983) rather than that of the star. This scenario has gained support over the last few years with the detection of HCOOCH₃ and CH₃OCH₃ in cold low-mass cores away from the YSO where the molecules are either released by shocks (Arce et al. 2008; Öberg et al. 2010) or by photodesorption (Bacmann et al. 2012) although the precise formation mechanisms are still unclear. Gas-phase processes may also contribute in some cases. Again, high spatial resolution and sensitivity such as ALMA will be key to testing this scenario.

6. Summary and conclusions

We have carried out a partial submillimeter line survey, targeting several complex organic species toward selected massive YSOs with strong evidence of circumstellar accretion disks, IRAS 20126+4104, IRAS 18089-1732 and G31.41+0.31. This is the first time that molecular abundances are reported for these well-known sources. The analysis is performed using both the rotation diagram method and spectral modeling. The inferred rotation temperatures and molecular abundances are compared to sources without reported disk structures analyzed using the same techniques. The molecules can be divided into two different groups based on their rotation temperatures, independent of source type. In particular, the cold (<100 K) species have remarkably constant rotation temperatures from source to source. The warm (>100 K) species exhibit more scatter possibly due to optical depth effects, non-LTE conditions and radiative excitation. The column densities peak for the same sources, with G31.41+0.31 being chemically the richest of the studied sources.

The relative abundances of species follow the same pattern, and no chemical differentiation could be established between the two source types within the (considerable) uncertainties. The lack of chemical differentiation between massive YSOs with and without observed disk-like structure suggests similarity in the physical conditions and level of UV exposure for all sources, for example through outflow cavities. Alternatively, the complex molecules may already be formed in the cold prestellar stage under similar conditions. Ultimately and indirectly, these results imply that the mechanism of the formation of massive stars does not differ significantly from source to source. Future high angular and high sensitivity observations of optically thin lines will be key to distinguishing the different scenarios and derive accurate abundance ratios that can be directly compared with chemical models.

---

Acknowledgments

Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by a Spinoza grant from the Netherlands Organisation for Scientific Research (NWO), and by the European Community’s Seventh Framework Programme FP7/2007-2013 under grant agreement 238258 (LASSIE) and 291141 (CHEMPLAN).

References

Online material

Appendix A: Detected lines per species for all sources

The line assignment and detection is based on Gaussian fitting with the following criteria: (i) the fitted line position has to be within ±1 MHz of the catalog frequency; (ii) the FWHM is consistent with the those given in Table 1; and (iii) the peak intensity has at least a S/N = 3. The errors on the integrated intensities are computed as follows.

The integrated main beam temperatures, , were obtained by Gaussian fits to the lines (Eq. (A.1)). (A.1)with

where T₀ is the peak intensity and ΔV is the FWHM of the line. The error, , is calculated from Eq. (A.2). (A.2)with

and

where rms is the root mean square amplitude of the noise in the spectral bin δv, cal is the calibration uncertainty of the telescope, and σ’s are the statistical errors on T₀ and ΔV from the Gaussian fits.

The errors on the integrated intensities derived from Eq. (A.2) include all statistical errors from the Gaussian fit. (A.3)For undetected transitions the upper limits were determined as 3σ limits (Eq. (A.4)) using: (A.4)where 1.2 is the coefficient related to the calibration uncertainty of 20%.

CLASS was used to determine the Gaussian fits and the uncertainties in the individual parameters. The formal errors on the integrated intensities derived from Eq. (A.3) in some cases yield a S/N < 2.5. This is caused by: (i) a conservative estimate of the statistical error on the FWHM parameter in the Gaussian fitter of CLASS; (ii) all statistical errors are included into our error calculation. Considering the higher S/N on the peak intensity, a more traditional error estimate without the statistical errors from the Gaussian fit (Eq. (A.3)) would result in a S/N > 3 for the integrated intensity as well. All weak line fits were confirmed by visual inspection.

Appendix B: Rotation diagrams

RTD diagrams for H₂CO, CH₃OH, C₂H₅OH, HNCO, NH₂CHO, CH₃CN, C₂H₅CN, HCOOCH₃, CH₃OCH₃, CH₂CO and CH₃CCH. Lines from different frequency bands are corrected for differential beam dilution assuming a R_{T = 100  K} source size for the warm species and 14″ source size for the cold species, as indicated in the plots. Optically thin, unblended lines with S/N ≳ 2 are marked with filled circles. Lines with high S/N (≲2) are marked with diamonds, those included in the fit with filled diamonds. Optically thick lines and blended lines are marked with open squares and triangles, respectively. Upper limits are marked with arrows. Upper limits used to constrain the fit are marked with filled stars. The error bars are calculated using Eq. (A.2).

Appendix C: Weeds model parameters

Appendix D: Additional detections

Several lines of other species were found in the observed frequency ranges, particularly in the line-rich G31.41+0.31. Table D.1 lists additional detections with column densities obtained from the Weeds analysis with the best T_ex when this could be derived. Several transitions of acetone, CH₃COCH₃, are detected. In G31, the best agreement is obtained at 100 K, whereas for IRAS 18089, this is T_ex = 300 K. No CH₃COCH₃ was detected in IRAS 20126, and the tabulated value is an upper limit at 300 K. Figure D.1 shows the strongest acetone lines in G31 together with the Weeds model.

Glycoaldehyde, HCOCH₂OH, has previously been detected in G31 (Beltrán et al. 2009), has several lines in the covered ranges. All lines are however blended with other transitions, and only upper limit could therefore be derived reliably. The tabulated column density is for a temperature of 300 K, constrained by non-detection of lines with low E_up.

Two transitions of HC₃N are detected at 218.325 GHz (E_up131.0 K) (J = 24 → 23) and 354.697 GHz (E_up340.5 K) (J = 39 → 38). The Weeds model on the two lines, well separated in E_up, gives a beam averaged column density of 7 × 10¹³ cm^-2 and a temperature of 150 K for G31. The HC₃N emission cannot be matched with emission contained within a 2.0″ volume as the low-E_uptransition becomes optically thick. For IRAS 20126, the tabulated column density is at 100 K, constrained by the two transitions. Only the high E_up transition was covered for IRAS 18089 giving a column density of 3 × 10¹³ cm^-2 at 150 K and 6 × 10¹³ cm^-2 at 100 K.

Several CH₃NH₂ transitions are observed with a T_ex of 150 K. Blends with other species make the rotation temperature and column density inaccurate, however. Several unidentified lines coincide with NH₂CH₂CH₂OH, but only upper limits of 1.0× 10¹⁶ cm^-2 at 100-300 K could be derived reliably using Weeds.

Several weak lines of NH₂CN are detected, with a T_ex of 50 K derived from the Weeds model, but line blends make the rotation temperature and column density inaccurate. For CN, SO, ³⁴SO, O¹³CS, OC³⁴S, SO₂, ³³SO₂, C³⁴S, and CP not enough lines were observed to derive T_ex values, and tabulated column densities are assuming a temperature of 50 K.

One line of HCN, J = 4–3, was observed. HCN emission in G31 is composed of a broad (12 km s^-1) and a narrow (6 km s^-1) component (redshifted by 2.5 km s^-1) causing self-absorption. For IRAS 18089, the broad component has a width of 7 km s^-1 and narrow 4 km s^-1 (redshifted by 1 km s^-1). For IRAS 20126 the HCN emission has a very broad 20 km s^-1 component blueshifted by 2 km s^-1 and a narrow (5 km s^-1) component redshifted by 1 km s^-1. The tabulated column densities are those derived from the H¹³CN column density.

All Tables

All Figures