Chemical complexity in high-mass star formation: An observational and modeling case study of the AFGL 2591 VLA 3 hot core

We present a detailed observational and modeling study of the hot core VLA 3 in the high-mass star-forming region AFGL 2591, which is a target region of the NOrthern Extended Millimeter Array (NOEMA) large program CORE. Using NOEMA observations at 1.37 mm with an angular resolution of ~0."42 (1 400 au at 3.33 kpc), we derived the physical and chemical structure of the source. We modeled the observed molecular abundances with the chemical evolution code MUSCLE (MUlti Stage ChemicaL codE). Results. With the kinetic temperature tracers CH3CN and H2CO we observe a temperature distribution with a power-law index of q = 0.41+-0.08. Using the visibilities of the continuum emission we derive a density structure with a power-law index of p = 1.7+-0.1. The hot core spectra reveal high molecular abundances and a rich diversity in complex molecules. The majority of the molecules have an asymmetric spatial distribution around the forming protostar(s), which indicates a complex physical structure on scales<1 400 au. Using MUSCLE, we are able to explain the observed molecular abundance of 10 out of 14 modeled species at an estimated hot core chemical age of ~21 100 years. In contrast to the observational analysis, our chemical modeling predicts a lower density power-law index of p<1.4. Reasons for this discrepancy are discussed. Conclusions. Combining high spatial resolution observations with detailed chemical modeling allows us to derive a concise picture of the physical and chemical structure of the famous AFGL 2591 hot core. The next steps are to conduct a similar analysis for the whole CORE sample, and then use this analysis to constrain the chemical diversity in high-mass star formation to a much greater depth.


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 > 8 M 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 Send offprint requests to: C. Gieser, e-mail: gieser@mpia.de 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 Hii (UC Hii) 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 nitrogenand 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 A v ≈ 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 C 70 . 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 (  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 in the lower left corner. The gray dashed circle indicates the primary beam of our NOEMA observations (23 ). 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 H 2 CO and CH 3 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, CH 3 OCHO, and CH 3 OCH 3 , 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. 2014Gerner et al. , 2015Urquhart 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. 2007c;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 ). In total, 18 HMS-FRs 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 H 2 O masers with the 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-infrared 2 , a large-scale outflow 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 Hii regions (Trinidad et al. 2003). The major source in the region is the VLA 3 hot core with a systemic velocity of v LSR = −5.5 km s −1 and a luminosity of 2×10 5 L suggesting a protostar with ∼40 M . 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 H 2 O 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, H 18 2 O, and SO 2 , Wang et al. (2012) found a Hubble law-like expansion within the inner 2 400 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 Hii 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 H 2 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 3 000 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 (H 2 S, 13 CS), double-peaked distributions (HC 3 N, OCS, SO, SO 2 ), and ring-like structures (CH 3 OH).
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 physicalchemical 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.

Observations
This study is part of the NOEMA large program CORE . 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 GHz 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. (subm.).
The following sources were observed for calibration during the AFGL 2591 tracks: the quasars 2013+370 and 2037+511 for phase and amplitude, a strong quasar (3C454.3) for bandpass, and the star MWC349 for absolute flux calibration. The data were calibrated with the CLIC package in gildas 3 . 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 deconvolved with the SDI algorithm (Steer et al. 1984), but see Mottram et al., subm., 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.

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 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 (1 400 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 H 2 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).

Line identification
An average NOEMA spectrum (spatially averaged over 1. 6×1. 6, indicated by the white rectangle in Fig. 2) around the continuum peak is used to identify emission lines detected at S /N > 5 (σ = 0.34 mJy beam −1 in the average spectrum). The molecular transitions are taken from the Cologne Database for Molecular Spectroscopy 4 (CDMS, Müller et al. 2005) and Jet Propulsion Laboratory database 5 (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 in the Appendix. We detect many transitions from simple molecules (e.g., SO, HNCO, OCS), but also a forest of COM lines, among which CH 3 OH and CH 3 CN have the strongest emission lines. For some molecules (e.g., HC 3 N and CH 3 OH) 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 A&A proofs: manuscript no. 35865corr   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 compact emission 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).

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.

Temperature structure
The molecules H 2 CO and CH 3 CN 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 CH 3 CN emission is strong around hot cores and can trace dense hightemperature gas (Fuente et al. 2014), whereas H 2 CO is used to trace the temperature of the colder extended envelope (Mangum & Wootten 1993). The upper state energies E u /k B of the observed transitions of H 2 CO and CH 3 CN are listed in Table A.1.
To derive the physical parameters of the molecular gas such as column density N and rotation temperature T rot , we fit the spectral lines with XCLASS 6 (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 H 2 density of ∼10 7 cm −3 in the central part within a radius of ∼2 500 au. The critical density is < 10 7 cm −3 for most of the observed molecules (CH 3 OH 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 T ex = T rot = T kin . 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 temperature T rot , column density N, linewidth ∆v, and velocity offset from the systemic velocity v off ) 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 model parameter, 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 CH 3 CN and H 2 CO emission lines from the merged data with one emission component in XCLASS, we derive the temperature structure assuming T rot = T kin adopting a threshold of 10σ above the noise level (σ = 1.5 mJy beam −1 = 0.23 K). For H 2 CO, we fit the 3 0,3 − 2 0,3 and 3 2,1 − 2 2,1 lines. For CH 3 CN, we fit the K = 4, 5, and 6 components of the J = 12 − 11 ladder, including the weak K = 0 − 3 of CH 13 3 CN transitions adopting an isotopic ratio of 60 (Wilson & Rood 1994). We exclude the K = 0 − 3 components of CH 3 CN 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 CH 3 CN, an algorithm chain using the Genetic (300 iterations) and the Levenberg-Marquardt (50 iterations) algorithms worked best. For the H 2 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 factor 7 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 H 2 CO and CH 3 CN temperature maps are shown in Fig. 4(a). The temperature distribution of CH 3 CN 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 H 2 CO shows more extended emission tracing the colder envelope with lower temperatures ranging 7 η = θ 2 source θ 2 source +θ 2 beam from 30 K to 150 K. In the case of H 2 CO, we restrict the analysis to the central area of ∼3 ×3 as the outer edges of the H 2 CO maps are sensitivity limited. For CH 3 CN 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 H 2 CO temperature distribution reaches a noisy plateau at temperatures > 100 K. This is due to the fact that the observed H 2 CO 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 K and 68 K), the fitting parameters N and T rot 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. 4(b) and there is a consistent overlap between the derived temperatures of both molecules at ∼2 000 au. The inner part below 700 au is spatially unresolved in our observations. The slope of the H 2 CO temperature profile gets steeper in the outer region at radii > 4 000 au. These data points are not included in the fit, as a large fraction of the H 2 CO line fluxes are below the adopted threshold of 10σ with a clear deviation from spherical symmetry as seen in Fig. 4(a). By combining and fitting (minimum χ 2 method) the CH 3 CN and H 2 CO temperature profiles, we derive the radial temperature profile from radii ranging from 700 au to 4 000 au: Palau et al. (2014) determined the temperature structure using a completely independent method fitting simultaneously the spectral energy distribution (SED) and radial intensity profiles.  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 H 2 CO temperature map. (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 H 2 CO (green data points) and CH 3 CN (blue data points). A power-law fit was performed on spatial scales in the range 700 − 4 000 au to determine the 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 temperature profile 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.
These authors find a value of q = 0.40 ± 0.01 and a temperature of 250 ± 20 K at 1 000 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.

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 H 2 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. 2007b). 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) 10 1 10 2 uv Distance [k ] 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 ∼3 000−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 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., Beuther et al. 2002;Mueller et al. 2002;Hatchell & van der Tak 2003;Beuther et al. 2007b;Zhang et al. 2009) and a detailed discussion is given in Sect. 5.1.

H 2 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(H 2 ) can be determined (Hildebrand 1983): Here S ν is the flux density, m H 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 cm 2 g −1 (for dust grains with a thin icy mantle at 1.3 mm at a gas density of 10 6 cm −3 ; Table 1 in Ossenkopf & Henning 1994). In each pixel, T is taken from the CH 3 CN and H 2 CO temperature maps shown in Fig. 4 The total gas mass can be derived as where F ν is the integrated flux density and d denotes the distance (3.33 ± 0.11 kpc, Rygl et al. 2012).
The resulting H 2 column density distribution based on Eq. 4 is shown in Fig. 6. The H 2 column density reaches 1.5×10 24 cm −2 in the center. Assuming a spherical source with a diameter of ∼5 000 au (1.5 ) along the line of sight as traced by our data in Fig. 2, the average H 2 volume density in the center is 2×10 7 cm −3 .
Using Eq. 5, the integrated total gas mass of the hot core core is 6.9±0.5 M 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 ), 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 × 10 4 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.

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 H 2 O 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 Array 8 (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 6 5 − 5 4 line obtained by the merged (NOEMA + IRAM 30 m) data tracing scales down to ∼1 400 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.   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 H 2 O maser bow-shocks are only 14 years, 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.

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 line data, is used to quantify the line emission properties of all detected molecules: rotation temperature T rot , column density N, linewidth ∆v, and velocity offset v off from the systemic velocity. 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 S /N of the HC 3 N;v 5 =1/v 7 =3, CH 2 CO, NH 2 CHO, CH 3 COCH 3 , and C 2 H 3 CN 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 HC 13 CCN and 13 CH 3 OH 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 simulta- neously in a single fit. Then the rotation temperature, linewidth, and velocity offset are modeled together, and the column density is calculated according to the assumed isotopic ratios. The following isotopologues were fitted simultaneously with their main species: O 13 CS, 34 SO 2 , H 13 2 CO, CH 13 3 CN. We apply fixed isotopic ratios taken from Wilson & Rood (1994): 12 C/ 13 C = 60 and 32 S/ 34 S = 22. As these authors do not discuss the 32 S/ 33 S isotopic ratio, the 33 S 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 v off from the systemic velocity (v LSR = −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 in the Appendix. 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 13 CO). 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 CH 3 CN and HNCO). The median temperature uncertainty is ∼14 % and the median column density uncertainty is ∼22 %. We calculate the average 32 S/ 33 S isotopic ratio using SO 2 and 33 SO 2 to be 211 ± 58. Using CS emission, Chin et al. (1996) found that for an ensemble of star-forming regions 32 S/ 33 S = 153 ± 40, which is in agreement with our results. The 32 S/ 33 S isotopic ratio in the solar system has a lower value of 127 (Anders & Grevesse 1989). Table 3. Positions of the peak integrated intensity taken from the data shown in Fig. 8 and Fig. 9 and classification of the morphology.  tion. 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 lineblending is not an issue. These maps clearly indicate chemical differentiation: a classification of the spatial distribution is sum-marized in Table 3 peak, sulfur-and nitrogen-bearing species are most abundant where the density and temperature are highest (e.g., SO, OCS, SO 2 , HC 3 N, CH 3 CN). Emission from most COMs is strongest in a ring-like structure towards the northeast, most visible in CH 3 OH. The emission of CH 2 CO 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 (O 13 CS, H 13 2 CO, HC 13 CCN) 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 H 13 2 CO, CH 3 OH, CH 2 CO, and CH 3 OCH 3 are consistent with a shock model coupled to a large gas-grain chemical network for a timescale of ∼2 000 years. 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.

Spatial molecular distribution
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.

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(Wakelam et al. , 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 r out have a constant inner part for r < r in and obey a power-law profile for r ≥ r in : The parameters of the physical structure (r in , r out , T 0 , p) can either be fixed or be varied as fit parameters within given ranges.

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 100 mag 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 . 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 CH 3 OH (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

Parameter
Value or Range Radiation field: CR ionization rate ζ CR 5(−17) s −1 X-ray ionization rate ζ X 0 s −1 Extinction at the cloud edge A v 100 mag UV photodesorption yield 1(−05) Grain properties: Grain radius r g 0.1 µm Dust density ρ d 3 g cm −3 Gas-to-dust mass ratio η 150 Surface diffusivity E Diff /E Bind 0.4 Mantle composition Mg x Fe x−1 SiO 4 (olivine) Pre-HMC stages IRDC stage a Inner radius r in 12 700 au Outer radius r out 0.5 pc Temperature T 0 at r in 11.3 K Temperature power-law index q 0.0 (isothermal) Density n 0 at r in 1.4(5) cm −3 Density power-law index p 1.5 τ IRDC 16 500 years HMC stage: Inner radius r in 50 − 700 au Inner radius step ∆r in 50 au Outer radius r out 7 000 − 31 000 au Outer radius step ∆r out 2 000 au Temperature power-law index q 0.41 Temperature T  an amorphous (unstructured) silicate composition. Each dust particle has 1.88×10 6 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 H 2 and CO is included as well.
For consistency with our observational analysis, we use a gas-todust 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 years 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 7 000 au (Fig. 2). From the observed radial temperature profile shown in Fig. 4(b), 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 in Sect. 3.3.1 (Eq. 1) and thus for each physical model the temperature T 0 at the inner radius is calculated according to T 0 = 255 K × ( r in 691 au ) −0.41 . In total, we run the chemical code ALCHEMIC on 2912 physical models (#r in × #r out × #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 years on 40 radial grid points.

Best-fit model
Using a χ 2 analysis, the observed column densities N obs are compared to the computed model column densities N mod 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 τ > 10 4 years, 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 magnitude 9 . Table 5. Best-fit model of the AFGL 2591 VLA 3 hot core.
Best-fit Model Parameters Inner radius r in 700 au Outer radius r out 29 000 au Temperature power-law index q 0.41 Temperature T 0 at r in 253.7 K Density power-law index p 1.0 Density n 0 at r in 1.7 (7)  The column densities derived from the observations given as input are listed in Table 2. The H 2 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 9 In the current MUSCLE version, uncertainties of the observed column densities cannot yet be included as an additional constraint. from the observed C 18 O abundance assuming an isotopic ratio of 16 O/ 18 O = 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 ≈ 4 600 years, 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 N model N obs = 40, 48, and 18, respectively, while for H 2 CO the modeled column density is too low with N model N obs = 3.1 × 10 −3 . Palau et al. (2017) found that there is an enhancement of gas-phase H 2 CO sputtered off the grains along an outflow of the high-mass protostar IRAS 20126+4104. Thus, one explanation for the underestimation of the H 2 CO abundance in our model could be that our chemical model does not include shock chemistry. The integrated intensity of H 13 2 CO (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 (r in × r out × p)-parameter space is shown in Fig. 11 for all models with χ 2 < 0.56. The bestfit model lies in a parameter space region with p ≈ 1.0 − 1.4 and τ HMC ≈ 3 000 − 6 000 years. With MUSCLE we cannot constrain r in or r out 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 H 2 CO abundances stay approximately constant. The sulfur-bearing species SO, OCS, and SO 2 have a high abundance in the center up to ∼4 000 au (> 124 K), with a sharp drop afterwards. The nitrogen-bearing species HNCO, HC 3 N, and CH 3 CN are enhanced in the central region up to 3 000 au (> 140 K) as well. A similar profile is seen for t-HCOOH, CH 3 OH, and CH 3 OCHO, but the abundances drop at 6 000 au (> 100 K). In contrast, the abundance of C 2 H 5 CN rises only towards the outer region > 10 000 au (> 85 K). For C 2 H 5 CN, higher model tem- peratures would be required in order to form it efficiently in the densest inner region. With a hot core stage timescale of τ HMC ≈ 4 600 years, the total lifetime (τ = τ IRDC + τ HMC ) of AFGL 2591 VLA 3 is 21 100 years. This is in rough agreement with previous studies that determined the chemical age of the hot core τ of a few 10 4 years: ∼30 000 years (Doty et al. 2002), ∼50 000 years (Stäuber et al. 2005), ∼80 000 years (Doty et al. 2006), and 10 000 − 50 000 years (Kaźmierczak-Barthel et al. 2015). Theoretical models of the formation of massive stars also predict lifetimes < 10 5 years for luminous YSOs Mottram et al. 2011).

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 H 2 CO and CH 3 CN and derived the temperature power-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 nonhomologous 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 observations of 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. (subm.) 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 C 34 S 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 largescale envelope may be flatter.

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 forming protostar. 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 A&A proofs: manuscript no. 35865corr Table 6. Literature comparison of the density power-law index p.
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 age would be > 49 000 years 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 1 100 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 form during 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 to study 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. (2014Gerner et al. ( , 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(Trinidad et al. , 2013Sanna 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 Fig. 8 and Fig. 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.

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 physicalchemical 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.37 mm 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 13 C, D, 18 O, 33 S, and 34 S and observed pure rotational (v=0), vibrational (up to v=3), and torsional (v t =1) transitions.
-The molecular thermometers CH 3 CN and H 2 CO were used to derive the temperature structure. We obtain a power-law index of q = 0.41 ± 0.08. While the observed CH 3 CN lines (E u /k B = 183, 247, and 326 K) trace hot and compact emission, the distribution of the observed H 2 CO emission lines (E u /k B = 21 and 68 K) are more extended tracing the lowtemperature gas. By combining both temperature profiles, we were able to determine the kinetic temperature over a wide range of radii from 700 to 4 000 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 H 2 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.5 M . 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 < 1 400 au or that we might be observing a disk wind. -Using XCLASS the column density N and rotation temperature T rot 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 CH 3 CN trace the central high-temperature region. Integrated intensity maps of detected molecular lines reveal that sulfurand nitrogen-bearing molecules have the strongest emission towards the 1.37 mm continuum peak (e.g., SO, OCS, HC 3 N, CH 3 CN, NH 2 CHO, C 2 H 3 CN, C 2 H 5 CN). Complex molecules such as CH 3 OH, CH 3 OCHO, and CH 3 COCH 3 peak towards the north and northeast with an asymmetric ring-like structure. The integrated intensity of less abundant isotopologues (O 13 CS, H 13 2 CO, HC 13 CCN) 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 p model = 1.0. The hot core chemical age is τ ≈ 21 100 years 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.