Open Access
Issue
A&A
Volume 631, November 2019
Article Number A142
Number of page(s) 21
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201935865
Published online 08 November 2019

© C. Gieser et al. 2019

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Open Access funding provided by Max Planck Society.

1 Introduction

The formation of the most massive stars is an active field of research (for reviews, see, e.g., Beuther et al. 2007a; Bonnell 2007; Zinnecker & Yorke 2007; Smith et al. 2009; Tan et al. 2014; Schilke 2015; Motte et al. 2018). Since high-mass star-forming regions (HMSFRs) are less abundant than their low-mass counterparts and are typically located at large distances (> 1 kpc), it is challenging to study individual collapsing objects.

The majority of high-mass star formation (HMSF) takes place in the densest regions of molecular clouds. Based on observational and theoretical considerations, HMSF can be divided into several evolutionary stages (Beuther et al. 2007a; Zinnecker & Yorke 2007): the formation of massive stars begins in infrared dark clouds (IRDCs) harboring potentially short-lived high-mass starless cores and low- to intermediate-mass protostars (e.g., Pillai et al. 2006; Rathborne et al. 2006; Sanhueza et al. 2012; Zhang et al. 2015) and proceeds to form high-mass protostellar objects (HMPOs) with M > 8M showing gas accretion and molecular outflows (e.g., Beuther et al. 2002; Motte et al. 2007). During the hot molecular core (HMC) stage, also called the hot core stage, the central protostar(s) warm the surrounding dense envelope to temperatures higher than 100− 300 K within the innermost several thousand au and rich molecular spectra at (sub)mm wavelengths are observed (e.g., Belloche et al. 2013; Sánchez-Monge et al. 2017; Beltrán et al. 2018). At later stages, ultra-compact H II (UC H II) regions form where the protostars ionize their surrounding envelope, which is observed in strong free–free emission at cm wavelengths (e.g., Garay & Lizano 1999; Palau et al. 2007; Qin et al. 2008; Klaassen et al. 2018).

Numerous millimeter and submillimeter wavelength observations towards hot cores reveal a rich chemistry in nitrogen- and oxygen-bearing molecules (for a review, see Herbst & van Dishoeck 2009) as well as sulfur-bearing species (Herpin et al. 2009). In addition, complex organic molecules (COMs) are very abundant in those environments. Following the definition by Herbst & van Dishoeck (2009), a COM consists of six or more atoms. Molecules can form in the interstellar medium (ISM), where the gas and dust are shielded from strong stellar radiation (at a visual extinction of Av ≈ 200 mag in the dense parts of HMSFRs; Tan et al. 2014). A summary of the detected molecules in the ISM is given in McGuire (2018) ranging from diatomic species such as CO up to C70. The chemical composition of these objects can be modeled using a chemical kinetics approach (Semenov et al. 2010; Cuppen et al. 2017)and an appropriate network of astrophysically relevant chemical reactions, for example, the KInetic Database for Astrochemistry (KIDA, Wakelam et al. 2015) and the UMIST Database for Astrochemistry (UDfA, McElroy et al. 2013). The reaction data are obtained by quantum chemical calculations and ultra-high vacuum laboratory experiments (for a review, see Cuppen et al. 2017).

In the early phases of star formation at low temperatures of ~10 K, atoms and small molecules (e.g., CO) produced in the gas-phase stick to the grain surface creating icy mantles (Tielens & Hagen 1982; Hasegawa et al. 1992). Hydrogenation of CO leads to species such as H2 CO and CH3 OH on the grains. Garrod & Herbst (2006) have shown that reactions between radicals on dust grains can form the most widely detected complex molecules in hot cores, such as HCOOH, CH3 OCHO, and CH3 OCH3, more efficiently during a gradual warm-up phase compared to the formation in the gas-phase. Jets and molecular outflows launched from protostars and interacting with the surrounding material create shocks that temporarily increase the temperature and density (Draine & McKee 1993; Bachiller 1996; Zhang et al. 2005; McKee & Ostriker 2007; Anderl et al. 2013). In such shocked regions molecular dissociation (but also reactions with activation barriers) can occur, and molecules formed on the grain mantles can desorb into the gas-phase (e.g., an increased abundance of COMs can be observed; Palau et al. 2017). Silicon is mostly locked on grains, but due to powerful outflows it is sputtered off the grains and forms SiO in the gas-phase (Herbst et al. 1989; Schilke et al. 1997; Gerin 2013).

The molecular richness of HMSFRs makes them ideal astrochemical laboratories. However, comprehensive chemical surveys of a large sample of high-mass star-forming cores with a high angular resolution targeting the molecular gas still need to be carried out and studied (e.g., Beuther et al. 2009; Wang et al. 2014; Feng et al. 2016). Surveys using single-dish telescopes have been made (e.g., Jackson et al. 2013; Gerner et al. 2014, 2015; Urquhart et al. 2019), but their angular resolution is not sufficient to resolve individual cores and the resulting molecular abundances are beam-averaged over a large region, which may contain signals of chemically distinct cores at different evolutionary stages. Most chemical studies at high angular resolution focus on individual cores (e.g., Beuther et al. 2007b; Fuente et al. 2014; Beltrán et al. 2018), whereas the influence of the environment on the chemical evolution may play an important role. For example, Feng et al. (2016), Allen et al. (2017), and Mills et al. (2018) observed and investigated chemical differentiation between adjacent hot cores.

In the framework of a sample of HMSFRs observed with the same setup and high spatial resolution, we carried out the NOrthern Extended Millimeter Array (NOEMA)1 large program CORE at 1.37 mm (Beuther et al. 2018). In total, 18 HMSFRs were observed with an angular resolution of ~0″. 4. The broad frequency bandwidth of ~4 GHz covers a large number of molecular lines and thus we are able to study the chemical properties of the objects in this sample. A detailed description of the observations is given in Sect. 2. As a test case for the chemical survey, we start with a case study of AFGL 2591 that will be complemented in the future with an analysis of the other sources in the CORE sample (Gieser et al., in prep.).

The high-mass star-forming region AFGL 2591, located at α(J2000) = 20:29:24.8 and δ(J2000) = +40:11:19.6, has been thoroughly studied in the past. By measuring the trigonometric parallax of H2 O masers withthe Very Long Baseline Array (VLBA), Rygl et al. (2012) determined a distance of 3.33 ± 0.11 kpc. A multiwavelength overview of the region is shown in Fig. 1. In the near-infrared2, a large-scaleoutflow lobe can be observed in the east–west direction (also Fig. 13 in Hodapp et al. 2003). At centimeter wavelengths a cluster of at least five sources is revealed with the Very Large Array (VLA, Campbell 1984; Trinidad et al. 2003; Johnston et al. 2013). VLA 1 and VLA 2 are H II regions (Trinidad et al. 2003). The major source in the region is the VLA 3 hot core with a systemic velocity of vLSR = −5.5 km s−1 and a luminosity of 2 × 105 L suggesting a protostar with ~40 M (Sanna et al. 2012). VLA 4 is a ~9 M star with spectral type B2, while the faint object VLA 5 coincides with the position of the infrared source MASS 20292393+4011105 (Johnston et al. 2013). Traced by H2O maser emission, Sanna et al. (2012) and Trinidad et al. (2013) found another protostellar object, VLA 3-N, very close to the VLA 3 hot core (~0″. 4 north).

Our analysis is focused on the hot core AFGL 2591 VLA 3 whose kinematic and chemical properties have been studied at multiple wavelengths in the past. Based on the kinematic properties using PdBI observations of HDO, HO, and SO2, Wang et al. (2012) found a Hubble law-like expansion within the inner 2400 au, and argue that this is caused by a disk wind, a phenomenon that has been recently observed for G17.64+0.16 (Maud et al. 2018). The centimeter continuum emission was studied by van der Tak & Menten (2005) with the VLA. These authors modeled the emission assuming that the expansion of the H II region is halted by accretion flows from the molecular envelope, which reproduces the observed structure of AFGL 2591 VLA 3. Infrared observations with the Stratospheric Observatory for Infrared Astronomy (SOFIA) reveal molecular emission in H2 O and CS from the inner hot dense region (Indriolo et al. 2015; Barr et al. 2018). Jiménez-Serra et al. (2012) reported chemical segregation on scales smaller than 3000 au due to the interplay between molecular dissociation by UV radiation, ice evaporation, and high-temperature gas-phase chemistry. These authors identify three types of spatial distribution: molecules with a compact distribution at the location of the continuum peak (H2S, 13CS), double-peaked distributions (HC3N, OCS, SO, SO2), and ring-like structures (CH3OH).

We use AFGL 2591 VLA 3 as a case study to investigate the chemical complexity in high-mass star-forming regions based on observations at high angular resolution and using a physical-chemical model. We infer the physical structure of the source and compare it to theoretical predictions. Using the rich molecular line information, we not only quantify the chemical content, but also study the spatial distribution of the emission to address the following questions: which molecular emission traces the compact core? Which species show extended large-scale emission? Which molecules have an asymmetric distribution? Which of those molecules could be chemically linked? What is the modeled chemical age of the hot core?

This paper is organized as follows. The CORE observations are described in Sect. 2 with a summary of the data calibration and properties of the final data products. The observational results are reported in Sect. 3, where the physical and chemical structure are analyzed. The applied physical-chemical model of the source is presented in Sect. 4. The observational and chemical modeling results are discussed in Sect. 5. In addition, the results are discussed in the context of previous studies of the target region. The main conclusions and an outlook are given in Sect. 6.

thumbnail Fig. 1

Overview of the AFGL 2591 star-forming region. The background color map shows an infrared image obtained with Gemini in the K′ band. The green contours show the 3.6 cm radio continuum emission taken from Johnston et al. (2013). The contours are drawn at levels of 3, 5, 10, 20, 40, 80, and 160σ (σ = 0.03 mJy beam−1). The VLA beam size (0″. 43 × 0″. 40) is shown inthe lower left corner. The gray dashed circle indicates the primary beam of our NOEMA observations (23″).

Open with DEXTER
Table 1

Data products of AFGL 2591 used in this study.

2 Observations

This study is part of the NOEMA large program CORE (Beuther et al. 2018). The main goals of the project are to study of the fragmentation and disk formation processes, as well as the outflows and chemistry in high-mass star formation. A detailed description of the full survey, technical setup, scientific goals, and analysis of the continuum data is given in Beuther et al. (2018). The observations of the sources in the 1 mm band were carried out with NOEMA from 2014−2017 in the A, B, and D configurations.

Using the wide-band correlator WideX, spectra from 217 to 221 GHz with a spectral resolution of 1.95 MHz (~2.7 km s−1 at 1 mm) were obtained in H and V linear polarizations. In order to study the kinematics, high spectral resolution units were placed around the wavelengths of key molecular transitions (see Ahmadi et al. 2018; Beuther et al. 2018). A study of the kinematics of the CORE regions, including AFGL 2591, will be given in Ahmadi et al. (in prep.).

In order to include short-spacing information, single-dish observations with the IRAM 30 m telescope were also performed. The IRAM 30 m observations were carried out in four 4 GHz broad spectral windows centered at ~215, 219, 231, and 234 GHz with a spectral resolution of 0.2 MHz (~0.3 km s−1 at 219 GHz) in H and V linear polarizations with a beam size of ~12″. The calibration of the single-dish data and merging process of the interferometric and single-dish data are described in detail in Mottram et al. (2018).

The following sources were observed for calibration during the AFGL 2591 tracks: the quasars 2013+370 and 2037+511 for phase and amplitude, astrong quasar (3C454.3) for bandpass, and the star MWC349 for absolute flux calibration. The data were calibrated with the CLIC package in gildas3. The 1.37 mm continuum was extracted from line-free channels in the spectral line data and then subtracted from the spectral line data. The NOEMA data (line and continuum) were deconvolved with the Clark CLEAN algorithm (Clark 1980) as the extended emission is filtered out. The NOEMA data merged with the IRAM 30 m observations (hereafter merged data) were deconvolvedwith the SDI algorithm (Steer et al. 1984), but see Mottram et al. (2018), for a comparison of these two deconvolution methods. The NOEMA and merged data were smoothed to a common spectral resolution of 3 km s−1. The properties of the data products are shown in Table 1.

3 Observational analysis

3.1 Continuum

The 1.37 mm continuum emission of the AFGL 2591 VLA 3 hot core is shown in Fig. 2. Other VLA sources (VLA 1, VLA 2, and VLA 4) are not detected within the NOEMA primary beam at a root mean square (rms) sensitivity of 0.59 mJy beam−1. The central emission of VLA 3 has an approximately spherically symmetric distribution around the position of the peak flux. With an angular resolution of 0″. 4 (1400 au at 3.33 kpc) we resolve the envelope, but no disk structure. Around the core there is an elongation towards the northwest, and there are several globules towards the northeast and southeast with emission > 5σ, both of which could potentially be associated with the molecular outflow (as studied by, e.g., Gueth et al. 2003). The main contribution to the 1.37 mm continuum emission is dust emission; based on VLA observations of H2 O masers in the AFGL 2591 star-forming region, at wavelengths ≥1 cm the emission is affected by free–free emission from the ionized gas (Trinidad et al. 2003). At frequencies ≥ 43 GHz the emission is dominated by the dust (van der Tak & Menten 2005).

thumbnail Fig. 2

AFGL2591 continuum emission at 1.37 mm. The green contours show the 1.37 mm continuum emission with levels at −3 (dashed), 3, 5, 10, 20, 40, and 80σ (σ = 0.59 mJy beam−1). The beam size (0″. 46 × 0″. 36) is shown inthe lower left corner. The green star indicates the position of the 1.37 mm continuum peak. The white rectangle shows the region within which an average spectrum was obtained in the spectral line data for line identification (Sect. 3.2).

Open with DEXTER

3.2 Line identification

An average NOEMA spectrum (spatially averaged over 1″. 6 × 1″. 6, indicated bythe white rectangle in Fig. 2) around the continuum peak is used to identify emission lines detected at SN > 5 (σ = 0.34 mJy beam−1 in the average spectrum). The molecular transitions are taken from the Cologne Database for Molecular Spectroscopy4 (CDMS, Müller et al. 2005) and Jet Propulsion Laboratory database5 (JPL, Pickett et al. 1998). The spectrum with labeled molecular lines is shown in Fig. 3. The properties of a total of 67 detected lines are summarized in Table A.1. We detect many transitions from simple molecules (e.g., SO, HNCO, OCS), but also a forest of COM lines, among which CH3OH and CH3CN have the strongest emission lines. For some molecules (e.g., HC3N and CH3OH) in addition to the commonly detected rotational states, rotational transitions from vibrationally or torsionally excited states are also present.

Multiple spectral line surveys of the AFGL 2591 region have been carried out in the past covering partially 80− 360 GHz using the James Clerk Maxwell Telescope (JCMT) and the IRAM 30 m telescope (Bisschop et al. 2007; van der Wiel et al. 2011); and 480−1900 GHz with the Herschel space observatory (Ceccarelli et al. 2010; van der Wiel et al. 2013; Kaźmierczak-Barthel et al. 2014). Bisschop et al. (2007) and Ceccarelli et al. (2010) reported that AFGL 2591 has a “line-poor” spectrum compared to other massive protostars, but with beam sizes >10″ the compactemission around the hot core cannot be resolved sufficiently and is heavily beam diluted. Our 1.37 mm spectral line data at subarcsecond resolution shows that the hot core is indeed a common hot core with a line-rich spectrum (Fig. 3).

3.3 Physical environment of AFGL 2591

We use the high angular resolution dust continuum (NOEMA) and merged (NOEMA + IRAM 30 m) spectral line data of AFGL 2591 at 1.37 mm to analyze the physical structure of the source.

3.3.1 Temperature structure

The molecules H2CO and CH3CN can probe the kinetic temperature when multiple transitions at different energy levels are observed (e.g., Zhang et al. 1998; Rodón et al. 2012). The CH3CN emission is strong around hot cores and can trace dense high-temperature gas (Fuente et al. 2014), whereas H2 CO is used to trace the temperature of the colder extended envelope (Mangum & Wootten 1993). The upper state energies EukB of the observed transitions of H2CO and CH3CN are listed in Table A.1.

To derive the physical parameters of the molecular gas such as column density N and rotation temperature Trot, we fit the spectral lines with XCLASS6 (eXtended Casa Line Analysis Software Suite, Möller et al. 2017). The XCLASS software models molecular lines by solving the 1D radiative transfer equation assuming local thermal equilibrium (LTE) conditions and an isothermal source. In dense sources, such as hot cores, the LTE assumption is a reasonable approximation for rotational lines; in Sect. 3.3.3 we estimate an average H2 density of ~107 cm−3 in the central part within a radius of ~2500 au. The critical density is < 107 cm−3 for most of the observed molecules (CH3OH being the only exception; see Table A.1). Therefore, the population of the molecules can be described by a Boltzmann distribution with a single characteristic temperature Tex = Trot = Tkin. In XCLASS, the line profiles are assumed to be Gaussian, but optical depth effects are included. Each molecule can be described by multiple emission and absorption components. The transition properties (e.g., Einstein coefficients, state degeneracies, partition functions) are taken from an embedded SQlite3 database containing entries from CDMS and JPL using the Virtual Atomic and Molecular Data Centre (VAMDC, Endres et al. 2016).

The XCLASS package provides new functions for CASA to model interferometric and single-dish data. By using the myXCLASSFit function, the user can fit the model parameters (source size θsource, rotation temperatureTrot, column density N, linewidth Δv, and velocity offset from the systemicvelocity voff) to observational data by using different optimization algorithms. In contrast, the myXCLASSMapFit function can be used to fit one or more complete data cubes. After fitting the selected pixel spectra, the myXCLASSMapFit function creates FITS images for each fitted modelparameter, where each pixel corresponds to the value of the optimized parameter for this pixel. The modeling can be done simultaneously with corresponding isotopologues and rotational transitions from vibrationally excited states. The ratio with respect to the main species or ground state can be either fixed or be introduced as an additional fit parameter.

Modeling the CH3CN and H2CO emission lines from the merged data with one emission component in XCLASS, we derive the temperature structure assuming Trot = Tkin adopting a threshold of 10σ above the noise level (σ = 1.5 mJy beam−1 = 0.23 K). For H2 CO, we fit the 30,3−20,3 and 32,1 −22,1 lines. For CH3CN, we fit the K = 4, 5, and 6 components of the J = 12 − 11 ladder, including the weak K = 0−3 of CHCN transitions adopting an isotopic ratio of 60 (Wilson & Rood 1994). We exclude the K = 0−3 components of CH3CN in order to reduce errors due to the high optical depth of those transitions, which may cause self-absorption (further discussed in Appendix B and in Ahmadi et al. 2018).

A detailed description of the optimization algorithms in XCLASS is given in Möller et al. (2013). It is necessary to carefully check whether the applied algorithms converge towards local minima or towards the aimed global minimum. For CH3CN, an algorithm chain using the Genetic (300 iterations) and the Levenberg–Marquardt (50 iterations) algorithms worked best. For the H2 CO fitting only the Levenberg-Marquardt algorithm (50 iterations) was used. In order to reduce the number of fit parameters from five to four, the source size θsource is fixed to 4″ for both molecules. In this case the beam filling factor7 is ~1 and therefore the emission is assumed to be resolved with the NOEMA beam θbeam ≈ 0″. 42. We show in Sect. 3.5 that the detected molecular emission is resolved in our data.

The derived H2CO and CH3CN temperature maps are shown in Fig. 4a. The temperature distribution of CH3CN is approximately spherically symmetric and is centered around the 1.37 mm continuum peak tracing the dense inner core with temperatures up to ~200 K, while H2CO shows more extended emission tracing the colder envelope with lower temperatures ranging from 30 to 150 K. In the case of H2 CO, we restrict the analysis to the central area of ~3″ × 3″ as the outer edges of the H2CO maps are sensitivity limited. For CH3CN this is not an issue as the emission is less extended. A discussion of the uncertainties of the fit parameters is given in Sect. 3.4. In the central region, the H2CO temperature distribution reaches a noisy plateau at temperatures > 100 K. This is due to the fact that the observed H2CO lines cannot trace the hottest gas (as seen in Fig. 7 in Rodón et al. 2012). As the rotation temperature is determined from the line ratios, its estimate would be more accurate if more than three transitions at different energy levels could have been observed. With only two lines at different upper energy levels (21 and 68 K), the fitting parameters N and Trot are degenerate.

Assuming spherical symmetry, we derive the radial temperature profile of both maps with the position of the 1.37 mm continuum peak emission as the center at α(J2000) = 20:29:24.88 and δ(J2000) = +40:11:19.47. The radial temperature profiles are shown in Fig. 4b and there is a consistent overlap between the derived temperatures of both molecules at ~2000 au. The inner part below 700 au is spatially unresolved in our observations. The slope of the H2 CO temperature profile gets steeper in the outer region at radii >4000 au. These data points are not included in the fit, as a large fraction of the H2 CO line fluxes are below the adopted threshold of 10σ with a clear deviation from spherical symmetry as seen in Fig. 4a. By combining and fitting (minimum χ2 method) the CH3CN and H2CO temperatureprofiles, we derive the radial temperature profile from radii ranging from 700 to 4000 au: (1)

Palau et al. (2014) determined the temperature structure using a completely independent method fitting simultaneously the spectral energy distribution (SED) and radial intensity profiles. These authors find a value of q = 0.40 ± 0.01 and a temperature of 250 ± 20 K at 1000 au for AFGL 2591 (their Tables 1 and 2). In a recent infrared study by Barr et al. (2018), hot abundant CS is detected with a temperature of 714 ± 59 K from a region <130 au. This is consistent with the temperature profile from our study, where this temperature is reached at a radius of 56 au.

thumbnail Fig. 3

Average spectrum (spatially averaged over 1″. 6 × 1″. 6, indicated by the white rectangle in Fig. 2) of the AFGL 2591 hot core labeled with molecular emission lines with SN > 5 (indicated by the red horizontal line, σ = 0.34 mJy beam−1 in the average spectrum). Unidentified lines are labeled UL. A detailed table of the properties of the molecular transitions is given in Appendix A.

Open with DEXTER
thumbnail Fig. 4

Panel a: H2CO and CH3CN temperature maps derived with XCLASS at a threshold of 10σ. The black contours show the 1.37 mm continuum emission with levels at 5, 10, 20, 40, and 80σ (σ = 0.59 mJy beam−1). The gray dash-dotted circle indicates the outer radius of the radial temperature fit (bottom panel). The beam size (0″. 47 × 0″. 36) is shown in the lower right corner of the H2CO temperaturemap. Panel b: radial temperature profiles obtained from the data in the top panel. The data are binned in steps of half a beam size (~700 au) for H2 CO (green data points) and CH3CN (blue data points). A power-law fit was performed on spatial scales in the range 700− 4000 au to determinethe temperature power-law index q. The gray dash-dotted line indicates the outer radius of the radial temperature fit. At radii smaller than 700 au the temperatureprofile cannot be determined due to the finite beam size (gray shaded area). The transparent data points were excluded from the fitting. The best fit is shown by the red solid line. The error of the fit (one standard deviation) is shown by the red dashed lines.

Open with DEXTER
thumbnail Fig. 5

Vector-averaged complex visibilities V of the 1.37 mm continuum as a function of uv distance. A power-law fit (red solid line) was performed in order to derive the power-law index α. The gray data points were excluded from the fit. The error on the fit (one standard deviation) is shown by the red dashed lines.

Open with DEXTER

3.3.2 Density structure

The peak flux density of the 1.37 mm continuum is 55.47 mJy beam−1 corresponding to a brightness temperature of 8.5 K (Fig. 2), whereas the gas has temperatures > 200 K (Fig. 4). We can therefore assume that the 1.37 mm continuum emission is optically thin and can thus be used as a tracer of the H2 column density which we analyze in the next section. However, because the continuum data are missing short-spacings, the images cannot be used to derive reliable intensity profiles. To overcome this problem, the intensity profiles can be fit in the Fourier-plane (e.g., Beuther et al. 2007c). The density profile can be determined using the previously derived temperature power-law index q and the measured continuum intensity distribution derived from the visibilities V (s) as a function of the uv distance s. Following the theoretical work of Adams (1991) and Looney et al. (2003), the density power-law index p can be derived assuming spherical symmetry: (2)

The vector-averaged complex visibilities of the 1.37 mm continuum emission, computed using the miriad software (Sault et al. 1995), are plotted against the uv distance in Fig. 5. The visibilities were shifted such that the position of the 1.37 mm continuum peak is at the phase center at α(J2000) = 20:29:24.88 and δ(J2000) = +40:11:19.47. With a power-law fit (minimum χ2 method), an index of α =−0.87 ± 0.05 is obtained on a physical scale from ~3000−67 000 au. Visibilities with a uv distance s > 200 kλ are excluded from the fit due to an increase in the noise level. Combined with the previous determination of the temperature index q (Sect. 3.3.1), the density index p can be calculated according to Eq. (2) as (3)

We obtain a density power-law index of p = 1.7 ± 0.1. This result agrees with Palau et al. (2014), where these authors find p = 1.80 ± 0.03 for AFGL 2591. Observations of HMSFRs find similar values for the density power-law index (e.g., van der Tak et al. 2000; Beuther et al. 2002, 2007c; Mueller et al. 2002; Hatchell & van der Tak 2003; Zhang et al. 2009) and a detailed discussion is given in Sect. 5.1.

thumbnail Fig. 6

Column density distribution of H2 derived fromthe 1.37 mm continuum emission at a 5σ threshold. The green contours show the 1.37 mm continuum emission with levels at 5, 10, 20, 40, and 80σ (σ = 0.59 mJy beam−1). The beam size (0″. 46 × 0″. 36) is shown in the lower left corner.

Open with DEXTER

3.3.3 H2 column density and gas mass

Optically thin dust emission dominates in the millimeter and submillimeter regime as discussed in the previous section. As noted in Sect. 3.1, the contribution from the ionized gas around AFGL 2591 can be neglected at 1.37 mm. Assuming that the gas and dust are thermally coupled, the molecular hydrogen column density N(H2) can be determined (Hildebrand 1983): (4)

Here Sν is the flux density, mH the mass of a hydrogen atom, Ω the beam solid angle, and Bν (T) the Planck function. Throughout the paper, we assume a gas-to-dust mass ratio of η = 150 (Draine 2011), a mean molecular weight of μ = 2.8, and a dust opacity of κν = 0.9 cm2 g−1 (for dust grains with a thin icy mantle at 1.3 mm at a gas density of 106 cm−3; Table 1 in Ossenkopf & Henning 1994). In each pixel, T is taken from the CH3CN and H2CO temperature maps shown in Fig. 4a with T = T(CH3CN) if T(CH3CN) > T(H2CO) and T = T(H2CO) elsewhere. Due to the temperature and density gradients, this is a rough approximation as CH3CN traces the inner dense region, while H2CO traces the colder envelope.

The total gas mass can be derived as (5)

where Fν is the integrated flux density and d denotes the distance (3.33 ± 0.11 kpc, Rygl et al. 2012).

The resulting H2 column density distribution based on Eq. (4) is shown in Fig. 6. The H2 column density reaches 1.5 × 1024 cm−2 in the center. Assuming a spherical source with a diameter of ~5000 au (1.5″) along the line of sight as traced by our data in Fig. 2, the average H2 volume density in the center is 2 × 107 cm−3.

Using Eq. (5), the integrated total gas mass of the hot core core is 6.9 ± 0.5M for emission >5σ within the innermost 10 000 au, assuming that the error budget is dominated by the uncertainties of the distance and flux density. A missing flux of 84% due to spatial filtering (Beuther et al. 2018), means that this should be taken as a lower limit. With continuum SED modeling of submillimeter JCMT SCUBA observations, van der Tak et al. (2013) estimate a gas mass of 373 M within a radius of 7.1 × 104 au. Based on SCUBA observations with a 22″. 9 beam at 850 μm (Di Francesco et al. 2008), the gas mass of the complete large-scale clump is estimated to be 638 M (using the same values for κν and η, and an average temperature of the region of T = 69 K derived from the IRAM 30 m data) within an effective radius of 77″. 6 corresponding to 1.3 pc (Table 1 in Beuther et al. 2018). This suggests that the AFGL 2591 hot core is embedded in a large-scale envelope providing a large gas reservoir.

3.3.4 Outflow(s)

Infrared observations reveal a large-scale east–west outflow from the VLA 3 hot core as shown in Fig. 1. In this section, we investigate how this compares to the outflow morphology traced by our high angular resolution observations. The outflow structure of AFGL 2591, as traced on various spatial scales in CO, SiO, and SO emission, is presented in Fig. 7. The properties of the different data products are summarized in Table 1.

The left panel of Fig. 7 shows the CO 2− 1 emission obtained by the IRAM 30 m telescope tracing scales down to ~37 000 au. There is a clear east–west(redshifted–blueshifted) bipolar outflow with the position of the 1.37 mm continuum peak in the center. The emission 30″ north originates from the large-scale filament in which the AFGL 2591 star-forming region is embedded. The observed outflow structure of AFGL 2591 VLA 3 is consistent with subarcsecond resolution observations of the same line by Jiménez-Serra et al. (2012), infrared observations (shown in Fig. 1), CO 3− 2 observations by Hasegawa & Mitchell (1995), and H2O maser observations by Sanna et al. (2012).

Another outflow tracer is SiO, which is produced by shocks (Schilke et al. 1997). No SiO line is present in the WideX frequency setup of the CORE observations with NOEMA. Archival Submillimeter Array8 (SMA) observations of the SiO 5−4 line at 217.105 GHz (project number 2010B-S108) are presented in the central panel in Fig. 7 tracing scales down to ~21 000 au (Jiménez-Serra, priv. comm.). The redshifted emission peaks at the position of the 1.37 mm continuum peak, while the blueshifted part peaks towards the north having an elongated structure towards the northwest direction. The integrated intensity around the central line shown by the gray scale reveals an extended structure towards the east. Within the SiO data at intermediate angular resolution, no clear outflow direction can be identified.

The right panel of Fig. 7 shows the outflow structure traced by the SO 65 − 54 line obtained by the merged (NOEMA + IRAM 30 m) data tracing scales down to ~1400 au. The redshifted outflow emission peaks towards the 1.37 mm continuum peak is similar to what is seen in SiO with a faint elongation towards the west. The blueshifted outflow part shows a clear direction towards the east in contrast to the large-scale outflow observed in CO where the blueshifted outflow points to the opposite direction.

The outflow has different morphologies using several outflow tracers at different spatial scales as shown in Fig. 7. While the large-scale east–west(redshifted–blueshifted) outflow traced by CO agrees with previous observations, at subarcsecond resolution the SO line might trace an east–west(blueshifted–redshifted) outflow. A similar outflow morphology is found for the intermediate-mass hot core IRAS 22198+6336 with two bipolar outflows: one of the outflows is traced by SiO emission, while the other outflow is traced by HCO+ emission (Sánchez-Monge 2011). While we do not resolve individual fragmented collapsing objects in the central region, different outflow directions traced by different molecules may indicate that AFGL 2591 VLA 3 hosts multiple chemically distinct young stellar objects (YSOs) with outflows. The outflow of the VLA 3-N object is traced by bow-shock structures of H2 O maser emission and has a north–south direction (Sanna et al. 2012; Trinidad et al. 2013). This source could be a good candidate to explain the anomalous direction of the SiO emission compared to the direction of the more evolved CO outflow. The inferred timescales for the H2 O maser bow-shocks are only 14 yr, as inferred by Sanna et al. (2012), which is consistent with the idea that this second outflow might be very young and hence we see it only in SiO. The emission might also be caused by a disk wind, as studied by Wang et al. (2012), for the AFGL 2591 hot core. Observations of outflow tracers at a higher angular resolution are required in order to investigate the outflow morphology in greater detail, which will be possible in the future by NOEMA.

thumbnail Fig. 7

Outflow structure of AFGL 2591 at different spatial scales traced by IRAM 30 m observations of CO (left panels), SMA observations of SiO (central panels), and NOEMA + IRAM 30 m observations of SO (right panels). Upper panels: background gray scale shows the integrated intensity of the central line emission. The blue and red contours show the blue- and redshifted line emission. The beam size is given in the lower left corner in each panel. The green star marks the position of the 1.37 mm continuum peak. The orange box indicates the change of field of view. Lower panels: spectra of the outflow tracers extracted from the position of the 1.37 mm continuum peak (green star). The systemic velocity of −5.5 km s−1 is indicated by the green dashed line. The blue and red dashed lines show the integration intervals for the blue- and redshifted line emission. The integration interval of the central line emission is set between the two ranges.

Open with DEXTER

3.4 Molecular content towards the continuum peak

The spectrum towards the position of the 1.37 mm continuum peak, extracted from the merged (NOEMA + IRAM 30 m) spectral linedata, is used to quantify the line emission properties of all detected molecules: rotation temperature Trot, column density N, linewidth Δv, and velocity offset voff from the systemicvelocity. A detailed description of the XCLASS fitting procedure and error estimation using myXCLASSFit is given in Appendix B. In contrast to the core-averaged spectrum shown in Fig. 3, the signal-to-noise ratio (S/N) of the HC3N;v5 = 1/v7 = 3, CH2CO, NH2CHO, CH3COCH3, and C2 H3CN emission is too low to fit the lines with XCLASS in the spectrum extracted at the position of the 1.37 mm continuum peak. The HC13CCN and 13CH3OH lines are blended and are therefore excluded from the fitting as well. We show in Sect. 3.5 that the molecular emission is spatially resolved, hence we adopt a beam filling factor η of 1 by fixing the source size θsource to 4″. This also reduces the fit parameter set from five to four. With XCLASS, it is also possible to model isotopologues simultaneously in a single fit. Then the rotation temperature, linewidth, and velocity offset are modeled together, and the column density is calculated according tothe assumed isotopic ratios. The following isotopologues were fitted simultaneously with their main species: O13CS, 34SO2, HCO, CHCN. We apply fixed isotopic ratios taken from Wilson & Rood (1994): 12C/13C = 60 and 32S/34S = 22. As these authors do not discuss the 32S/33S isotopic ratio, the 33S isotopologues are fitted individually.

The results for the rotation temperature and column densities are given in Table 2. The spectral resolution of 3 km s−1 of the merged data is too low to study in detail the kinematic properties of the emission lines, but is sufficient to determine the rotation temperature and column density. The mean linewidth Δv of the detected molecules is 5.2 ± 0.6 km s−1 and the mean velocity offset voff from the systemic velocity (vLSR = −5.5 km s−1) is 0.0 ± 0.3 km s−1. A comparison between the observed and modeled best-fit spectrum is shown in Fig. B.1. If only one strong transition is present in the spectrum, there is a large degeneracy between the column density and rotation temperature, thus the uncertainties of the fit parameters are high (e.g., for 13CO). Due to the large optical depth of the SO line (bottom panel in Fig. B.1), the derived column density can only be taken as a lower limit. However, if many transitions with different energy levels of a molecule are detected, the modeling with XCLASS is more reliable (e.g., for CH3CN and HNCO). The median temperature uncertainty is ~14% and the median column density uncertainty is ~22%. We calculate the average 32S/33S isotopic ratio using SO2 and 33SO2 to be 211 ± 58. Using CS emission, Chin et al. (1996) found that for an ensemble of star-forming regions 32S/33S = 153 ± 40, which is in agreement with our results. The 32S/33S isotopic ratio in the solar system has a lower value of 127 (Anders & Grevesse 1989).

Table 2

XCLASS fitting results towards the 1.37 mm continuum peak position.

3.5 Spatial molecular distribution

We investigate the spatial distribution of the molecular emission with line integrated intensity maps at high angular resolution. Figures 8 and 9 shows the integrated intensities obtained from the NOEMA spectral line data. The flux is integrated along five channels covering a width of 15 km s−1 (the average line FWHM is 5 km s−1, Sect 3.4) centered at the rest frequency of each transition. We carefully chose transitions where line-blending is not an issue. These maps clearly indicate chemical differentiation: a classification of the spatial distribution is summarized in Table 3. As we observe a similar chemical segregation as reported by Jiménez-Serra et al. (2012), we adopt a similar classification for the morphology of the spatial distribution of the molecular emission. Type I molecules have a single-peaked distribution (spatially), while type II molecules have a double-peaked distribution (spatially). Type III molecules have a ring-like morphology. Towards the 1.37 mm continuum emission peak, sulfur- and nitrogen-bearing species are most abundant where the density and temperature are highest (e.g., SO, OCS, SO2, HC3N, CH3CN). Emission from most COMs is strongest in a ring-like structure towards the northeast, most visible in CH3OH. The emission of CH2CO is shifted ~0″. 4 north towards the position of VLA 3-N which might be in a colder and younger stage. At this position, the emission of COMs is also bright but more spread out. The asymmetric emission from less abundant and optically thin isotopologues (O13CS, HCO, HC13CCN) may trace shocked regions where their abundances are enhanced. Palau et al. (2017) show that such asymmetries in molecular emission were also found at similar spatial scales of ~600 au towards the prototypical high-mass protostar IRAS 20126+4104, and also that the observed asymmetries for HCO, CH3OH, CH2CO, and CH3OCH3 are consistent with a shock model coupled to a large gas-grain chemical network for a timescale of ~2000 yr. The integrated intensity of the unidentified line (UL) at 219.616 GHz (bottom right panel in Fig. 9) is centered around the position of the 1.37 mm continuum peak.

The complex emission of various types of species shows that the formation of molecules depends heavily on the local physical conditions, such as the envelope, shocked regions, and outflows. It is very likely that the source consists of multiple chemically distinct sub-components.

4 Chemical model

The observational analysis of the AFGL 2591 CORE data reveal a large degree of chemical complexity both in abundances and spatial distribution. In this section, we use the physical-chemical modeling code MUSCLE (MUlti Stage ChemicaL codE, a detailed description of the model is given in Gerner et al. 2014) to study the formation and destruction of the observed molecular abundances.MUSCLE computes a static radial physical structure coupled with the time-dependent chemical kinetics code ALCHEMIC (Semenov et al. 2010; Semenov 2017) including both gas-phase and grain-surface chemistry with reaction rates obtained from KIDA (Wakelam et al. 2012, 2015). The computed abundance profiles are beam-convolved and thus converted into model column densities as a function of time. The model column densities are compared to the observed column densities, and with a χ2 analysis the best-fit physical structure and chemical age are determined.

The temperature T(r) and the density profile n(r) of the spherically symmetric model core with outer radius rout have a constant inner part for r < rin and obey a power-law profile for rrin: (6)

The parameters of the physical structure (rin, rout, T0, p) can either be fixed or be varied as fit parameters within given ranges.

thumbnail Fig. 8

Integrated intensity maps of detected molecular emission lines. In each panel the molecule and transition are noted (see Table A.1). The color map shows the integrated intensity of the line emission. The white contours indicate the 5σ level of the integrated intensity (σ = 1.2 K km s−1). The green contours show the 1.37 mm continuum emission with levels at 5, 10, 20, 40, and 80σ (σ = 0.59 mJy beam−1). The beam size of the spectral line data (0″. 47 × 0″. 36) is shown in the lower left corner. The white bar in the top right corner indicates a spatial scale of 3330 au (1″).

Open with DEXTER

4.1 Initial abundances and conditions

The adopted MUSCLE input parameters are summarized in Table 4. We assume that internal sources do not generate X-rays: ζX = 0 s−1. High energetic UV and X-ray radiation emitted by the forming YSO(s) in the center can affect the envelope material of the AFGL 2591 hot core through the outflow cavity (as studied by, e.g., Benz et al. 2007), but in our simplistic spherically symmetric model, the high energetic X-rays would be immediately absorbed by the dense gas in the center. A high extinction of 100mag in the outer envelope shields the core from the interstellar radiation field. At infrared wavelengths, the AFGL 2591 region seems to be isolated within 30″, so external heating is not important (van der Tak et al. 2000). Only cosmic rays (CRs) can penetrate into the central part of the core. We employ a cosmic ray ionization rate of 5 × 10−17 s−1 (van der Tak & van Dishoeck 2000). We use a low value of 10−5 for the UV photodesorption yield, as measured for CH3OH (Cruz-Diaz et al. 2016; Bertin et al. 2016), which gives the probability that a CR-induced UV photon kicks a molecule off an icy grain.

Grain-surface reactions depend on the geometrical structure of dust particles (Gerin 2013). We consider spherical grains with an amorphous (unstructured) silicate composition. Each dust particle has 1.88 × 106 surface sites (Semenov et al. 2010) where species can freeze out and chemical reactions can occur via the Langmuir–Hinshelwood mechanism (Hasegawa et al. 1992; Biham et al. 2001). Neutral species and electrons can stick to the grains with a probability of 100%. Only hydrogen atoms are allowed to tunnel through the surface barriers. After a reaction occurs, products can directly desorb from the grain with a probability of 1% (Garrod et al. 2007; Vasyunin & Herbst 2013). In addition, species can desorb by thermal processes, CRs, and UV radiation. Self-shielding of H2 and CO is included as well. For consistency with our observational analysis, we use a gas-to-dust mass ratio of 150 (Draine 2011). The diffusivity of species on dust grains can be parameterized by the ratio of the diffusion energy to the binding energy and as it is poorly constrained, we employ an experimentally motivated value of 0.4 (Cuppen et al. 2017).

The evolutionary past of the physical structure and chemical composition of AFGL 2591 is not known. As initial abundances, we use the best-fit IRDC abundances from the study by Gerner et al. (2015). For 19 well-studied IRDCs, column densities of 18 molecules including deuterated species were obtained using IRAM 30 m observations at 1 and 3 mm. The column densities of the molecules were averaged and modeled to create abundances of a template IRDC. These authors derive an IRDC chemical age of τIRDC = 16 500 yr using MUSCLE (Table A.4 and A.8 in Gerner et al. 2015).

Finally, the physical structure of the AFGL 2591 hot core has to be defined. Our NOEMA data are missing the short-spacings in the continuum yielding a missing flux of 84%, so we can only determine a lower limit for the outer radius of 7000 au (Fig. 2). From the observed radial temperature profile shown in Fig. 4b, we can adopt an upper limit for the inner radius of 700 au which is half a beam size. The density index p is varied from 1.0 to 2.5. The temperature profile is fixed to the radial profile derived from the observations presented inSect. 3.3.1 (Eq. (1)) and thus for each physical model the temperature T0 at the inner radius is calculated according to . In total, we run the chemical code ALCHEMIC on 2912 physical models (#rin × #rout × #p = 14 × 13 × 16). For each model, it takes ~75 s using 32 cores to compute the chemical evolution with ALCHEMIC up to 100 000 yr on 40 radial grid points.

thumbnail Fig. 9

Integrated intensity maps of detected molecular emission lines. In each panel the molecule and transition are noted (see Table A.1). The color map shows the integrated intensity of the line emission. The white contours indicate the 5σ level of the integrated intensity (σ = 1.2 K km s−1). The green contours show the 1.37 mm continuum emission with levels at 5, 10, 20, 40, and 80σ (σ = 0.59 mJy beam−1). The beam size of the spectral line data (0″. 47 × 0″. 36) is shown in the lower left corner. The white bar in the top right corner indicates a spatial scale of 3330 au (1″).

Open with DEXTER
Table 3

Positions of the peak integrated intensity taken from the data shown in Figs. 8 and 9 and classification of the morphology.

Table 4

Model input parameters.

4.2 Best-fit model

Using a χ2 analysis, the observed column densities Nobs are compared to the computed model column densities Nmod for each physical model and time step. Vasyunin et al. (2004) studied the error propagation of the abundances of chemical models due to the large uncertainties and approximations of rate coefficients. These authors estimated an error of 0.5−1.5 orders of magnitude depending on the size of the molecule. Wakelam et al. (2005) found that for older hot cores with τ > 104 yr, the uncertainties in the chemical model can be large for some species and therefore care should be taken when comparing observed and modeled abundances. We therefore assume that a molecule is modeled well if the modeled and observed column densities agree within one order of magnitude9.

The column densities derived from the observations given as input are listed in Table 2. The H2 column density is calculated from the 1.37 mm continuum emission at the peak position, as analyzed in Sect 3.3.3. The CO column density is inferred from the observed C18O abundance assuming an isotopic ratio of 16O/18O = 500 (Wilson & Rood 1994). For molecules whose lines from the vibrational ground state and from vibrationally excited states were observed, an average value for the column density is given as input.

The best-fit parameters from the model with the lowest χ2 are summarized in Table 5. The best-fit model has a hot core lifetime τHMC ≈ 4600 yr, with a density power-law index of p = 1.0. A comparison between the modeled and observed column densities is shown in Fig. 10. For 10 out of 14 molecules the modeled and observed column densities agree within the uncertainties. For SO, DCN, and HNCO the model column density is too high compared to the observations with , 48, and 18, respectively, while for H2CO the modeled column density is too low with . Palau et al. (2017) found that there is an enhancement of gas-phase H2 CO sputtered off the grains along an outflow of the high-mass protostar IRAS 20126+4104. Thus, one explanation for the underestimation of the H2 CO abundance in our model could be that our chemical model does not include shock chemistry. The integrated intensity of HCO (Fig. 8) shows an enhancement along a potentially blueshifted NE outflow part as traced by SO (Fig. 7).

The χ2 and τHMC distribution in the (rin × rout × p)-parameter space is shown in Fig. 11 for all models with χ2 < 0.56. The best-fit model lies in a parameter space region with p ≈ 1.0−1.4 and τHMC ≈ 3000−6000 yr. With MUSCLE we cannot constrain rin or rout as within our model parameter space all values reach a small χ2 (Fig. 11).

Within the best-fit model it is possible to analyze the radial variation of the modeled abundances,which is shown in Fig. 12 at the time of the best-fit chemical age τHMC. The modeled CO, DCN, and H2CO abundances stay approximately constant. The sulfur-bearing species SO, OCS, and SO2 have a high abundance in the centerup to ~4000 au (> 124 K), with a sharp drop afterwards. The nitrogen-bearing species HNCO, HC3N, and CH3CN are enhanced in the central region up to 3000 au (>140 K) as well. A similar profile is seen for t-HCOOH, CH3OH, and CH3OCHO, but the abundances drop at 6000 au (>100 K). In contrast, the abundance of C2 H5CN rises only towards the outer region >10 000 au (>85 K). For C2 H5CN, higher model temperatures would be required in order to form it efficiently in the densest inner region.

With a hot core stage timescale of τHMC ≈ 4600 yr, the total lifetime (τ = τIRDC + τHMC) of AFGL 2591 VLA 3 is 21 100 yr. This is in rough agreement with previous studies that determined the chemical age of the hot core τ of a few 104 yr: ~30 000 yr (Doty et al. 2002), ~50 000 yr (Stäuber et al. 2005), ~80 000 yr (Doty et al. 2006), and 10 000−50 000 yr (Kaźmierczak-Barthel et al. 2015). Theoretical models of the formation of massive stars also predict lifetimes < 105 yr for luminous YSOs (Davies et al. 2011; Mottram et al. 2011).

Table 5

Best-fit model of the AFGL 2591 VLA 3 hot core.

thumbnail Fig. 10

Comparison of the observed (red) and modeled (blue) column densities.

Open with DEXTER
thumbnail Fig. 11

MUSCLE fit parameter space. Left panel: physical models with χ2 < 0.56. Right panel: best-fit chemical age τHMC of the physical models with χ2 < 0.56.

Open with DEXTER
thumbnail Fig. 12

Radial abundance profiles of the MUSCLE best-fit model, relative to the radial H2 density profile in the bottom right panel.

Open with DEXTER
Table 6

Literature comparison of the density power-law index p.

5 Discussion

5.1 Observed and modeled density and temperature structure of the AFGL 2591 hot core

In Sect.3.3 we present the analysis of the physical environment of the hot core based on the CORE observational data. First, we inferred the kinetic temperature structure using H2 CO and CH3CN and derived the temperaturepower-law index q = 0.41 ± 0.08 assuming spherical symmetry (Fig. 4), where the temperature profile may not follow a single power-law. Furthermore, using the visibilities of the 1.37 mm continuum emission (Fig. 5), we inferred the density power-law index p = 1.7 ± 0.1 of the source. Again, spherical symmetry was assumed for the analysis. In reality, this assumption is most likely not valid for both the temperature and density profile due to the outflow and the potential fragmentation on even smaller scales, as discussed in Sects. 3.3 and 3.5. However, the assumption of spherical symmetry is sufficient to describe the surrounding envelope.

Deriving the density structure in high-mass star-forming regions has not been attempted frequently, but a few studies have been put forward in previous years. Table 6 summarizes the results for the density profiles p obtained in other studies. There are several theoretical models that explain the collapse of a gas cloud into a core. Numerical calculations of a collapsing spherically symmetric cloud by Larson (1969) have shown that without initial pressure gradients the free-fall collapse is non-homologous and the density profile approaches p = 2.0. Shu (1977) studied the collapse of an isothermal cloud that is in hydrostatic equilibrium. The initial density distribution follows a power-law with p = 2.0 and the accretion rate is constant during the collapse. McLaughlin & Pudritz (1997) refined this simple model, considering a non-isothermal gas distribution and a logotropic equation of state, which can best explain observationsof both low- and high-mass cores from small to cloud scales. In this case, an initial value of p = 1.0 holds and during the collapse the accretion rate increases. In both models the gas eventually free-falls onto the center with p = 1.5. As summarized in Table 6, observational determinations of the density profile p range from 1− 2. One explanation for such large variations could be that during the evolution in HMSF the temperature and density structure may change drastically, especially when considering different formation and accretion mechanisms that require non-spherical symmetry. Furthermore, depending on the tracer and the spatial resolution, different spatial scales and regions are traced, which might have different density profiles and asymmetries.

Regarding the density profile of the AFGL 2591 hot core, our model predicts p < 1.4, in contrast to our observational value of p = 1.7 ± 0.1. We cannot determine whether the density profile is dominated by non-thermal processes (p = 1), free-fall collapse (p = 1.5), or thermal processes (p = 2). The observed temperature profile presented in Fig. 4 is not fully spherically symmetric and does not obey a single power-law. As a consequence, the model temperature in the outer region is higher than observed. As q and p are correlated, a steeper temperature index results in a flatter density index. Mottram et al. (2018) and Beltrán et al. (2018) observe steeper values for hot cores: q = 0.7 and q = 0.8, respectively. The density profile of AFGL 2591 was derived in a previous study based on single-dish observations: van der Tak et al. (1999) modeled CS and C34S emission lines to infer a density power-law index of 1.25 ± 0.25. Our observations of AFGL 2591 VLA 3 at high angular resolution suggest that p > 1.6. These differences indicate that the inner core may have a steep density profile, while the density profile of the large-scale envelope may be flatter.

5.2 Model uncertainties

The observed molecular column densities presented in Sect. 3.4 were applied to the physical-chemical 1D model MUSCLE in order to study the physical structure of the hot core and radial profiles of the molecular abundances. The observed molecular emission around the hot core (shown in Fig. 8 and Fig. 9) reveal that on subarcsecond scales the structure of the molecular gas is highly asymmetric. Therefore, the modeling approach of a spherically symmetric core is too simplistic to understand all the chemical properties on such small scales. For instance, dense gas may be shielded by the inner disk, molecules in the outflow cavity might be strongly affected by high-energy radiation of the formingprotostar. Unfortunately, a full 3D physical-chemical model is not yet available to investigate these effects further.

While MUSCLE is able to model several evolutionary stages of HMSF, within each stage the physical structure is static. There is no temperature and density evolution within a model stage, but in reality the core is collapsing, where the accretion disk and outflow evolve as well. The reaction rates are very sensitive to temperature and density changes and the high-temperature structure that we observe was preceded by a warm-up phase. We have no knowledge of the conditions of the pre-HMC stages of AFGL 2591. Therefore, the best-fit IRDC abundances from Gerner et al. (2015) were applied as a starting point for the HMC stage of AFGL 2591. We tested different initial conditions of the abundances based on previous MUSCLE studies: using the IRDC+HMPO stage from Gerner et al. (2015), the model performs worse with fewer molecules modeled well and a higher χ2 value. In that case the chemical agewould be >49 000 yr suggesting that the hot core is younger. Using the IRDC + HMPO stage from the CORE pilot region NGC 7538 S (Feng et al. 2016), we obtain similar results, but a steeper density profile. However, the modeled physical structure is much more compact with an outer radius of 1100 au and an extrapolation of the chemical abundances to our AFGL 2591 physical models with much higher outer radii is not reasonable. Many molecules, including simple COMs and their precursors readily formduring the cold IRDC phase (e.g., Vasyunina et al. 2014). Hence the initial atomic and molecular abundances play an important role. As studied by Wilson & Rood (1994) and Wannier (1980), among others, the abundances in the ISM depend on the local environment and there are gradients from the Galactic center towards the outer regions, as high-mass stars quickly enrich the local environment with heavy metals. Thus, it is necessary tostudy the chemical conditions in the ISM, the influence of the environment, and the evolutionary timescales in HMSFRs with a large statistical sample at a high spatial resolution, analogous to the work by Gerner et al. (2014, 2015), among others, who used single-dish observations. A study of the chemical content of all 18 regions at high angular resolution within the large program CORE will be presented in a future publication (Gieser et al., in prep.).

One uncertainty of the model is the assumption of the radiation field. During HMSF, when the star reaches the main sequence, it is still deeply embedded in the cloud and envelope. Therefore, the central source may already emit strong UV radiation and X-rays, which have a strong impact on the chemistry. This enhances the formation of ions but also the destruction of large molecules. The high-energy radiation can escape the central part of the source through the outflow and winds. Numerical simulations by Kuiper & Hosokawa (2018) have shown that photoionization feedback plays an important role in the outflows of massive stars on scales from 0.01−1 pc. Previous studies of AFGL 2591 have shown that some observed molecular abundances can only be explained when a central heating source is included in the model (Benz et al. 2007; Bruderer et al. 2009). The chemical segregation observed by Jiménez-Serra et al. (2012) was also explained by molecular UV photo-dissociation and high-temperature gas-phase chemistry in the central region. In addition, as recent observations by Csengeri et al. (2018) suggest, the physical structure of high-mass envelopes can be impacted by additional processes, such as shocks due to accretion.

Grain-surface reactions play an important role in the formation of molecules in the ISM (e.g., Potapov et al. 2017). In our model, we assumed spherical silicate dust grains with uniform size, but we neglect carbonaceous grains which may influence chemical reactions on the grains. Observations have revealed that dust grains follow a size distribution (Mathis et al. 1977; Weingartner & Draine 2001). Simpson et al. (2013) have shown that the dust grains around AFGL 2591 consist of both spherical and elongated grains. Both the grain size and form distribution define the effective surface area and, for larger areas, chemical reactions are more likely to occur.

In deriving the observed column densities and in the MUSCLE model, we assumed that AFGL 2591 VLA 3 consists of one hot core. Previous studies have claimed the presence of a second source (VLA 3-N) ~0″. 4 towards the north (Trinidad et al. 2003, 2013; Sanna et al. 2012). We observe several outflow directions in different outflow tracers as discussed in Sect. 3.3.4 and the chemical segregation of the molecular distribution shown in Figs. 8 and 9 indicates fragmentation on scales smaller than 1 400 au. Therefore, the observed molecular emission may come from several YSOs which we do not resolve. However, with our simplistic model we are still able to capture the most essential physics and chemistry of simple gas-phase species.

6 Conclusions

In this study, we analyzed continuum and spectral line data of AFGL 2591 VLA 3 at 1.37 mm obtained within the NOEMA large program CORE. We derived the radial temperature and density profiles, studied the outflow and determined the molecular abundances towards the position of the 1.37 mm continuum peak. We modeled the chemical abundances using the physical-chemical code MUSCLE and found a rough agreement with previous studies for the chemical age of the hot core. Our main conclusions of this study are the following:

  • The 1.37mm spectra reveal a large chemical complexity in O-, N-, and S-bearing species. In contrast to previous line surveys of the hot core at lower angular resolution, the spectra of the CORE data show that AFGL 2591 VLA 3 indeed has a typical hot core spectrum with many molecular emission lines. We found several isotopologues containing 13C, D, 18O, 33S, and 34S and observed pure rotational (v = 0), vibrational (up to v = 3), and torsional (vt = 1) transitions.

  • The molecular thermometers CH3CN and H2CO were used to derive the temperature structure. We obtain a power-law index of q = 0.41 ± 0.08. While the observed CH3CN lines (Eu/kB = 183, 247, and 326 K) trace hot and compact emission, the distribution of the observed H2CO emission lines (Eu/kB = 21 and 68 K) are more extended tracing the low-temperature gas. By combining both temperature profiles, we were able to determine the kinetic temperature over a wide range of radii from 700 to 4000 au.

  • The visibilities of the 1.37 mm continuum emission were used to derive the density power-law index. Assuming spherical symmetry and optically thin emission, we obtain a value of p = 1.7 ± 0.1. We derived the H2 column density from the 1.37 mm continuum emission and estimated the gas mass reservoir within the inner 10 000 au around the continuum peak to 6.9 ± 0.5M. This is a lower limit, as the NOEMA data miss a large amount of flux due to spatial filtering.

  • The emission of CO, SiO, and SO using single-dish and interferometric observations were used to study the outflow structure. Two potential blueshifted outflow directions can be identified, one towards the west and one towards the northeast. This suggests that there could be several fragmented cores on spatial scales <1400 au or that we might be observing a disk wind.

  • Using XCLASS the column density N and rotation temperature Trot of all detected molecules were determined. Some molecules lie in an outer, colder envelope such as DCN. This is also observed in another CORE target region NGC 7538 IRS9 (Feng et al., in prep.). Molecules such as CH3CN trace the central high-temperature region. Integrated intensity maps of detected molecular lines reveal that sulfur- and nitrogen-bearing molecules have the strongest emission towards the 1.37 mm continuum peak (e.g., SO, OCS, HC3N, CH3CN, NH2CHO, C2H3CN, C2H5CN). Complex molecules such as CH3OH, CH3OCHO, and CH3COCH3 peak towards the north and northeast with an asymmetric ring-like structure. The integrated intensity of less abundant isotopologues (O13CS, HCO, HC13CCN) may trace shocked regions. This chemical segregation can be explained by the fact that the formation of specific molecules is confined to different temperature and density regions (disk, outflow, envelope) with different dynamics (e.g., shocks), and that several fragmented cores may be hidden in the central part of the source.

  • We applied the physical-chemical model MUSCLE, including grain-surface and gas-phase chemistry, to the observed abundances of the hot core in order to derive the physical structure and chemical age. Our best-fit model can explain 10 out of 14 detected molecular species within the uncertainties with a density power-law index of pmodel = 1.0. The hot core chemical age is τ ≈ 21 100 yr consistent with previous studies. However, keeping the complicated structure of AFGL 2591 at scales < 10 000 au traced by the gas in mind, the assumption of spherical symmetry is too simplistic. With a static physical structure in each stage, molecules may not be modeled well due to the lack of a warm-up stage in MUSCLE.

The purpose of this work was to use the hot core AFGL 2591 VLA 3 as a case study in order to investigate the chemical abundance within the CORE spectral setup at 1.37 mm. The high spatial resolution and sensitivity obtained by NOEMA revealed a rich chemical complexity both in abundance and spatial distribution. The physical-chemical code MUSCLE is sufficient to model a large fraction of the observed molecular abundances obtained by interferometric observations at high angular resolution assuming spherical symmetry.

Acknowledgements

The authors would like to thank the anonymous referee whose comments helped improve the clarity of this paper. This work is based on observations carried out under project number L14AB with the IRAM NOEMA Interferometer and the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). C.G., H.B., A.A., and J.C.M. acknowledge support from the European Research Council under the Horizon 2020 Framework Programme via the ERC Consolidator Grant CSF-648505. D.S. acknowledges support from the Heidelberg Institute of Theoretical Studies for the project “Chemical kinetics models and visualization tools: Bridging biology and astronomy.” R.K. acknowledges financial support via the Emmy Noether Research Group on Accretion Flows and Feedback in Realistic Models of Massive Star Formation funded by the German Research Foundations (DFG) under grant no. KU 2849/3-1 and KU 2849/3-2. A.P. acknowledges financial support from UNAM-PAPIIT IN113119 grant, México. A.S.M. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG) via the Sonderforschungsbereich SFB 956 Conditions and Impact of Star Formation (subproject A6). I.J.-S. acknowledges partial support by the MINECO and FEDER funding under grants ESP2015-65597-C4-1 and ESP2017-86582-C4-1-R. This research made use of Astropy10, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

Appendix A Transition properties of detected molecular lines

Table A.1 shows the transition properties taken from the CDMS and JPL databases of all detected spectral lines in the average AFGL 2591 VLA 3 spectrum presented in Fig. 3. The critical density is calculated based on Einstein coefficients Aul and collisional rate coefficients Cul from the Leiden Atomic and Molecular Database11(LAMDA, Schöier et al. 2005): via the approximation (Shirley 2015).

Table A.1

Transition properties (rest frequency, species, quantum numbers, upper state energy, critical density) of detected molecular lines in the average AFGL 2591 spectrum shown in Fig. 3.

Vibrationally excited states are marked by vx (with x = 1, 2, 3,...). For CH3OH, torsional transitions (rotation of the methyl group within the molecule) are detected and marked with vt. If nothing is added next to the molecule name, the transition is pure rotational (v = 0). The quantum numbers are given depending on the symmetry of the molecule: JupperJlower for linear molecules, for symmetric top molecules, and for asymmetric top molecules.

Depending on the nuclear spin states of the three hydrogen atoms in the methyl group, CH3OH and CH3OCHO are denotedwith A (parallel proton spin) or E (anti-parallel proton spin), whereas CH3OCH3 and CH3COCH3 have two methyl groups (denoted with AA, EE, AE, or EA). Formic acid (HCOOH) exists in two conformers (c-HCOOH and t-HCOOH) depending on the orientation of the hydroxy (OH) group. Multiple transitions at the same frequency are blended within our spectral setup at a spectral resolution of 2.2 MHz.

Appendix B Spectral line modeling with XCLASS

Similar to the map fitting described in Sect. 3.3.1, we use XCLASS to fit the spectrum extracted towards the position of the 1.37 mmcontinuum peak. In order to find the best-fit parameters of all detected molecules with XCLASS, an algorithm chain with the Genetic (50 iterations) and Levenberg–Marquardt algorithm (50 iterations) is adopted (a detailed description of the algorithms is given in Möller et al. 2013). In the spectrum we estimate a rms noise of 0.26 K in a line-free frequency range at 217.664−217.706 GHz. For each molecule that is detected with emission lines SN > 5, the myXCLASSFit function models the non-blended emission lines with one emission component in order to find the best-fit parameters. In order to estimate the uncertainties of the fit parameters, we scale the intensity of spectrum to both + 20 and − 20%, which is the maximum expected uncertainty from the flux calibration of the CORE data, as discussed in Beuther et al. (2018). The two scaled spectra are likewise fitted with myXCLASSFit and the uncertainties are then estimated by the mean deviation from the best-fit parameters. A comparison between the observed and fitted spectrum is shown in Fig. B.1. There are large differences between the observed and modeled emission lines for molecules which have either non-Gaussian line profiles and/or a high optical depth (e.g., 13CO, SO, H2 CO and the CH3CN K = 3 transition).

Table B.1

Rotation temperatures calculated with RADEX.

The underlying assumption in the calculation of the synthetic spectra with XCLASS is that LTE conditions hold. In order to check that this assumption is valid, we use the non-LTE radiative transfer code RADEX12 (van der Tak et al. 2007) to calculate therotation temperature Trot of the molecules with the derived physical parameters as an input. We assume the CMB (cosmic microwave background) temperature for the background temperature (Tbg = 2.73 K). For the kinetic temperature, we use the CH3CN temperatureTkin = 200 K and column densities of each molecule derived towards the continuum peak. We assume the H2 volume density of n(H2) ≈ 107 cm−3 described in Sect. 3.3.3 and the linewidth of Δv = 5.2 km s−1 described inSect. 3.4. Local thermal equilibrium conditions are given when Trot = Tkin.

Table B.1 shows the computed rotation temperature for all molecular column densities towards the continuum peak that have available data in RADEX. For most molecules TexTkin, except for SO2 (9 K) and CH3OH (10 K). As vibrational excited SO2 with high upper energy levels (≳900 K) is also detected, part of the emission has to stem from a non-LTE environment (e.g., in a shocked region as studied by Jørgensen et al. 2004). The negative rotation temperatures of CH3OH denote that maser activity is expected and indeed hot cores are known to show strong CH3OH maser emission (e.g., Longmore et al. 2007).

thumbnail Fig. B.1

Comparison of the observed and fitted XCLASS spectrum. Upper panel: AFGL 2591 spectrum towards the continuum peak (black line) and the total XCLASS fit spectrum (red line). Middle panel: residuals (green line). Bottom panel: optical depth τ (blue line) of the transitions computed with XCLASS.

Open with DEXTER

References


1

Upgraded version of the Plateau de Bure Interferometer (PdBI).

2

Based on observations obtained at the Gemini Observatory acquired through the Gemini Observatory Archive, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), Ministério da Ciência, Tecnologia e Inovação (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).

7

.

8

The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica.

9

In the current MUSCLE version, uncertainties of the observed column densities cannot yet be included as an additional constraint.

All Tables

Table 1

Data products of AFGL 2591 used in this study.

Table 2

XCLASS fitting results towards the 1.37 mm continuum peak position.

Table 3

Positions of the peak integrated intensity taken from the data shown in Figs. 8 and 9 and classification of the morphology.

Table 4

Model input parameters.

Table 5

Best-fit model of the AFGL 2591 VLA 3 hot core.

Table 6

Literature comparison of the density power-law index p.

Table A.1

Transition properties (rest frequency, species, quantum numbers, upper state energy, critical density) of detected molecular lines in the average AFGL 2591 spectrum shown in Fig. 3.

Table B.1

Rotation temperatures calculated with RADEX.

All Figures

thumbnail Fig. 1

Overview of the AFGL 2591 star-forming region. The background color map shows an infrared image obtained with Gemini in the K′ band. The green contours show the 3.6 cm radio continuum emission taken from Johnston et al. (2013). The contours are drawn at levels of 3, 5, 10, 20, 40, 80, and 160σ (σ = 0.03 mJy beam−1). The VLA beam size (0″. 43 × 0″. 40) is shown inthe lower left corner. The gray dashed circle indicates the primary beam of our NOEMA observations (23″).

Open with DEXTER
In the text
thumbnail Fig. 2

AFGL2591 continuum emission at 1.37 mm. The green contours show the 1.37 mm continuum emission with levels at −3 (dashed), 3, 5, 10, 20, 40, and 80σ (σ = 0.59 mJy beam−1). The beam size (0″. 46 × 0″. 36) is shown inthe lower left corner. The green star indicates the position of the 1.37 mm continuum peak. The white rectangle shows the region within which an average spectrum was obtained in the spectral line data for line identification (Sect. 3.2).

Open with DEXTER
In the text
thumbnail Fig. 3

Average spectrum (spatially averaged over 1″. 6 × 1″. 6, indicated by the white rectangle in Fig. 2) of the AFGL 2591 hot core labeled with molecular emission lines with SN > 5 (indicated by the red horizontal line, σ = 0.34 mJy beam−1 in the average spectrum). Unidentified lines are labeled UL. A detailed table of the properties of the molecular transitions is given in Appendix A.

Open with DEXTER
In the text
thumbnail Fig. 4

Panel a: H2CO and CH3CN temperature maps derived with XCLASS at a threshold of 10σ. The black contours show the 1.37 mm continuum emission with levels at 5, 10, 20, 40, and 80σ (σ = 0.59 mJy beam−1). The gray dash-dotted circle indicates the outer radius of the radial temperature fit (bottom panel). The beam size (0″. 47 × 0″. 36) is shown in the lower right corner of the H2CO temperaturemap. Panel b: radial temperature profiles obtained from the data in the top panel. The data are binned in steps of half a beam size (~700 au) for H2 CO (green data points) and CH3CN (blue data points). A power-law fit was performed on spatial scales in the range 700− 4000 au to determinethe temperature power-law index q. The gray dash-dotted line indicates the outer radius of the radial temperature fit. At radii smaller than 700 au the temperatureprofile cannot be determined due to the finite beam size (gray shaded area). The transparent data points were excluded from the fitting. The best fit is shown by the red solid line. The error of the fit (one standard deviation) is shown by the red dashed lines.

Open with DEXTER
In the text
thumbnail Fig. 5

Vector-averaged complex visibilities V of the 1.37 mm continuum as a function of uv distance. A power-law fit (red solid line) was performed in order to derive the power-law index α. The gray data points were excluded from the fit. The error on the fit (one standard deviation) is shown by the red dashed lines.

Open with DEXTER
In the text
thumbnail Fig. 6

Column density distribution of H2 derived fromthe 1.37 mm continuum emission at a 5σ threshold. The green contours show the 1.37 mm continuum emission with levels at 5, 10, 20, 40, and 80σ (σ = 0.59 mJy beam−1). The beam size (0″. 46 × 0″. 36) is shown in the lower left corner.

Open with DEXTER
In the text
thumbnail Fig. 7

Outflow structure of AFGL 2591 at different spatial scales traced by IRAM 30 m observations of CO (left panels), SMA observations of SiO (central panels), and NOEMA + IRAM 30 m observations of SO (right panels). Upper panels: background gray scale shows the integrated intensity of the central line emission. The blue and red contours show the blue- and redshifted line emission. The beam size is given in the lower left corner in each panel. The green star marks the position of the 1.37 mm continuum peak. The orange box indicates the change of field of view. Lower panels: spectra of the outflow tracers extracted from the position of the 1.37 mm continuum peak (green star). The systemic velocity of −5.5 km s−1 is indicated by the green dashed line. The blue and red dashed lines show the integration intervals for the blue- and redshifted line emission. The integration interval of the central line emission is set between the two ranges.

Open with DEXTER
In the text
thumbnail Fig. 8

Integrated intensity maps of detected molecular emission lines. In each panel the molecule and transition are noted (see Table A.1). The color map shows the integrated intensity of the line emission. The white contours indicate the 5σ level of the integrated intensity (σ = 1.2 K km s−1). The green contours show the 1.37 mm continuum emission with levels at 5, 10, 20, 40, and 80σ (σ = 0.59 mJy beam−1). The beam size of the spectral line data (0″. 47 × 0″. 36) is shown in the lower left corner. The white bar in the top right corner indicates a spatial scale of 3330 au (1″).

Open with DEXTER
In the text
thumbnail Fig. 9

Integrated intensity maps of detected molecular emission lines. In each panel the molecule and transition are noted (see Table A.1). The color map shows the integrated intensity of the line emission. The white contours indicate the 5σ level of the integrated intensity (σ = 1.2 K km s−1). The green contours show the 1.37 mm continuum emission with levels at 5, 10, 20, 40, and 80σ (σ = 0.59 mJy beam−1). The beam size of the spectral line data (0″. 47 × 0″. 36) is shown in the lower left corner. The white bar in the top right corner indicates a spatial scale of 3330 au (1″).

Open with DEXTER
In the text
thumbnail Fig. 10

Comparison of the observed (red) and modeled (blue) column densities.

Open with DEXTER
In the text
thumbnail Fig. 11

MUSCLE fit parameter space. Left panel: physical models with χ2 < 0.56. Right panel: best-fit chemical age τHMC of the physical models with χ2 < 0.56.

Open with DEXTER
In the text
thumbnail Fig. 12

Radial abundance profiles of the MUSCLE best-fit model, relative to the radial H2 density profile in the bottom right panel.

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

Comparison of the observed and fitted XCLASS spectrum. Upper panel: AFGL 2591 spectrum towards the continuum peak (black line) and the total XCLASS fit spectrum (red line). Middle panel: residuals (green line). Bottom panel: optical depth τ (blue line) of the transitions computed with XCLASS.

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.