# Multiwavelength modelling of the circumstellar environment of the massive protostar AFGL 2591 VLA 3

F. A. Olguin, M. G. Hoare, K. G. Johnston, F. Motte, H. -R. V. Chen, H.

Beuther, J. C. Mottram, A. Ahmadi, C. Gieser, D. Semenov, et al.

► To cite this version:

F. A. Olguin, M. G. Hoare, K. G. Johnston, F. Motte, H. -R. V. Chen, et al. Multiwavelength modelling of the circumstellar environment of the massive protostar AFGL 2591 VLA 3. Monthly Notices of the Royal Astronomical Society, 2020, 498, pp.4721-4744. 10.1093/mnras/staa2406. insu-03636584

HAL Id: insu-03636584

https://insu.hal.science/insu-03636584

Submitted on 20 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Multiwavelength modelling of the circumstellar environment of the massive protostar AFGL 2591 VLA 3

F. A. Olguin<sup>1,2*</sup> M. G. Hoare<sup>2</sup> K. G. Johnston,<sup>2</sup> F. Motte,<sup>3,4</sup> H.-R. V. Chen,<sup>1</sup> H. Beuther,<sup>5</sup>

J. C. Mottram,<sup>5</sup> A. Ahmadi,<sup>5</sup> C. Gieser,<sup>5</sup> D. Semenov,<sup>5,6</sup> T. Peters,<sup>7</sup> A. Palau<sup>8</sup>, P. D. Klaassen<sup>9</sup>,

R. Kuiper,<sup>10</sup> Á. Sánchez-Monge<sup>11</sup> and Th. Henning<sup>5</sup>

Accepted 2020 August 5. Received 2020 August 5; in original form 2019 August 19

## ABSTRACT

We have studied the dust density, temperature, and velocity distributions of the archetypal massive young stellar object (MYSO) AFGL 2591. Given its high luminosity ($L = 2 \times 10^5 L_\odot$) and distance ($d = 3.3$ kpc), AFGL 2591 has one of the highest $\sqrt{L/d}$ ratio, giving better resolved dust emission than any other MYSO. As such, this paper provides a template on how to use resolved multiwavelength data and radiative transfer to obtain a well-constrained 2D axisymmetric analytic rotating infall model. We show for the first time that the resolved dust continuum emission from Herschel 70- µm observations is extended along the outflow direction, whose origin is explained in part from warm dust in the outflow cavity walls. However, the model can only explain the kinematic features from CH$_3$CN observations with unrealistically low stellar masses ($<15 M_\odot$), indicating that additional physical processes may be playing a role in slowing down the envelope rotation. As part of our three-step continuum and line fitting, we have identified model parameters that can be further constrained by specific observations. High-resolution mm visibilities were fitted to obtain the disc mass (6 $M_\odot$) and radius (2200 au). A combination of SED and near-infrared observations were used to estimate the luminosity and envelope mass together with the outflow cavity inclination and opening angles.

**Key words:** circumstellar matter – stars: formation – ISM: individual objects: AFGL 2591 – infrared: stars.

## 1 INTRODUCTION

High-mass stars play an important role in the evolution of galaxies, but their formation is still not well understood. Do high-mass stars form as a scaled-up version of isolated low-mass star formation (e.g. Krumholz, Klein & McKee 2012; Kuiper & Hosokawa 2018) or by some other means (e.g. competitive accretion, Smith, Longmore & Bonnell 2009; fragmentation-induced starvation, Peters et al. 2010; global hierarchical collapse, Vázquez-Semadeni et al. 2019)? To understand how high-mass stars form, we need to study the distribution of matter in star-forming regions, which is the result of the interaction of inflow and outflow processes driven by gravitational collapse, rotation, turbulence, magnetic fields, and radiation. The imprints left by these processes in the circumstellar material can be well studied in early stages of massive young stellar objects

(MYSOs), when the ionizing nature of the stellar radiation has not started to evaporate the core.

These MYSOs are radio weak (radio fluxes of few mJy, Hoare 2002), bright in the infrared (IR) and have luminosities $L \gtrsim 10^4 L_\odot$ (e.g. Mottram et al. 2011b). The radiation pressure on dust produced by such massive stars can be reduced by the high optical depths of an accretion disc (Kuiper et al. 2010), and also through bipolar outflow cavities opened by the interaction of jets and winds with the infalling material (e.g. Krumholz, McKee & Klein 2005; Cunningham et al. 2011; Kuiper, Yorke & Turner 2015; Kuiper, Turner & Yorke 2016). Temperatures lower than their zero-age main sequence counterparts are also thought to be responsible for the lack of ionizing radiation given their luminosity (Hoare & Franco 2007; Davies et al. 2011), which may be due to the forming star being swollen due to accretion (Hosokawa, Yorke & Omukai 2010; Kuiper & Yorke 2013; Palau et al. 2013; Haemmerlé & Peters 2016; Kuiper & Hosokawa 2018).

MYSOs are located within giant molecular clouds at larger distances (typical distances of 3 kpc) compared to nearby low-mass

<footnote n="1"><sup>1</sup>Institute of Astronomy and Department of Physics, National Tsing Hua University, Hsinchu 30013, Taiwan</footnote>

<footnote n="2"><sup>2</sup>School of Physics & Astronomy, E.C. Stoner Building, University of Leeds, Leeds LS2 9JT, UK</footnote>

<footnote n="3"><sup>3</sup>Institut de Planétologie et d'Astrophysique de Grenoble, Univ. Grenoble Alpes – CNRS-INSU, BP 53, F-38041 Grenoble Cedex 9, France</footnote>

<footnote n="4"><sup>4</sup>AIM Paris-Saclay, CEA/IRFU – CNRS/INSU – Univ. Paris Diderot, Service d'Astrophysique, CEA-Saclay, F-91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="5"><sup>5</sup>Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="6"><sup>6</sup>Department of Chemistry, Ludwig Maximilian University, Butenandtstr. 5-13, D-81377 Munich, Germany</footnote>

<footnote n="7"><sup>7</sup>Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany</footnote>

<footnote n="8"><sup>8</sup>Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, 58090 Morelia, México</footnote>

<footnote n="9"><sup>9</sup>UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="10"><sup>10</sup>Institute of Astronomy and Astrophysics, University of Tübingen, Auf der Morgenstelle 10, D-72076, Tübingen, Germany</footnote>

<footnote n="11"><sup>11</sup>I. Physikalisches Institut, Universität zu Köln, Zülpicher Str 77, D-50937 Köln, Germany</footnote>

<footnote>* E-mail: folguin@gapp.nthu.edu.tw</footnote>
star-forming clouds, and usually in clustered environments. Their study is therefore observationally challenging due to their location, the amount of gas and dust in their parental clouds and their fast formation ($10^4-3 \times 10^5$ yr time-scales; Mottram et al. 2011b; Russeil et al. 2010; Duarte-Cabral et al. 2013). Hence, to study these regions we must observe them at longer wavelengths and at high angular resolution. In the last decade, an effort to map large regions of the sky at high angular resolution has been undertaken using ground- (e.g. UKIRT,<fnref n="1" /> JCMT,<fnref n="2" /> APEX<fnref n="3" />) and space-based telescopes (e.g. Spitzer, Herschel) from near-IR to sub-millimetre wavelengths. These observations are also complemented by interferometric observations of specific sources at even higher resolutions, which allow the study of regions closer to the forming star (e.g. Maud et al. 2013; Beuther et al. 2018).

The study by Olguin et al. (2015) of Red MSX Source (RMS) survey (Lumsden et al. 2013) MYSOs as mapped by the Herschel IR Galactic Plane Survey (Hi-GAL; Molinari et al. 2010) shows that relatively isolated sources with high $\sqrt{L/d}$, with L is the source luminosity and d its distance, may be resolved at 70 µm. They also analysed three sources in the $l=30^\circ$ and $59^\circ$ fields by fitting their spectral energy distributions (SEDs) and 70-µm data with 1D radiative transfer models assuming a power-law density distribution. The density power-law index they obtained was shallower ($\le 1$) than expected for infalling material (index of 1.5). These results suggest that the far-IR emission may be dominated by warm dust from the outflow cavity walls rather than rotational flattening as suggested by earlier studies (e.g. de Wit et al. 2009), as the mapped emission is larger than the expected centrifugal radius. These findings are in line with the results of 2D axisymmetric radiative transfer modelling of mid-IR observations (e.g. de Wit et al. 2010).

In the light of these new data, a multiwavelength study of MYSOs can provide further constraints to the dust/gas density, temperature, and velocity distributions of their circumstellar matter. In this work, we present a method to fit multiwavelength data and assess the importance of specific observations to constrain these distributions. We apply this method to the well-studied MYSO AFGL 2591 to explain a large range of observations and compare the results with those in the literature.

## 1.1 The selected source: AFGL 2591

AFGL 2591 is a well-studied luminous ($L \sim 10^5 L_\odot$) source located at $3.3 \pm 0.1$ kpc (Rygl et al. 2012) in the Cygnus X sky region. Several radio continuum sources have been identified in this region: four sources have been classified as H II regions (VLA 1, 2, 4, 5; see Trinidad et al. 2003; Johnston et al. 2013), and one as a MYSO or young hot core (VLA 3; Trinidad et al. 2003; Gieser et al. 2019). An unknown source(s) associated with maser emission has also been detected (VLA 3-N; Trinidad et al. 2013). The MYSO (AFGL 2591 VLA 3) will be referred to as AFGL 2591 as it is the source that dominates the SED from the near-IR to millimetre wavelengths (Johnston et al. 2013). As such, this proto-typical MYSO/young hot core is one of the sources with the highest $\sqrt{L/d}$ in the RMS survey sample. It was covered by the Herschel/HOBYS survey<fnref n="4" /> (Motte et al. 2010). The HOBYS data for the Cygnus X region were published by

Schneider et al. (2016), however, there is no published close-up study of this source at 70 µm to date. The HOBYS survey was designed to map specific star-forming regions at lower scan speeds. Therefore, it is one of the best candidates to resolve the 70-µm emission since the observations are less subject to smearing effects compared to other Herschel observations such as Hi-GAL.

The presence of a jet within a large-scale outflow was inferred from early CO and HCO$^+$ molecular line observations (e.g. Lada et al. 1984; Hasegawa & Mitchell 1995), with the bipolar outflow cavity oriented in the E–W direction. The jet was also detected with radio interferometric observations at 3.6 cm, and has a position angle of ~100° (Johnston et al. 2013). The blue-shifted outflow cavity can also be observed in scattered light in the K band (2 µm, e.g. Preibisch et al. 2003), whilst shocked H$_2$ bipolar emission is detected in the same band (Poetzel, Mundt & Ray 1992; Tamura & Yamashita 1992). The blue-shifted cavity K-band emission presents several features (loops) formed probably by entrainment of material (e.g. Parkin et al. 2009). Hasegawa & Mitchell (1995) qualitatively constrained the full cavity opening angle to be <90° and its inclination with respect to the line of sight $\le 45^\circ$ from CO emission. Based on geometrical considerations from Gaussian fits to molecular line emission maps (e.g. SO$_2$), van der Tak et al. (2006) constrained the inclination angle to between $26^\circ$ and $38^\circ$. The near-IR polarization study by Minchin et al. (1991) found that a cone with an inclination of $35^\circ$ can reproduce their data, however, Simpson et al. (2013) found an inclination angle of $15^\circ$ by fitting the SED and then modifying this model in order to match the morphology of their HST polarization data. van der Tak et al. (1999) fitted a 1D spherically symmetric power-law density distribution to CS isotopologues observations and then modified their best model to include an empty bipolar cavity, and found that an opening angle of $60^\circ$ can better fit their data. However, the opening angle may be wider at the base of the cavity ($\theta_e \sim 100^\circ$) as implied by the spatial distribution and proper motion of water maser emission (Sanna et al. 2012). The detailed radiative transfer modelling by Johnston et al. (2013), whose objective was to fit the SED and JHK 2MASS slices along and perpendicular to the outflow direction, found that the inclination angle with respect to the line of sight is constrained to be between $30^\circ$ and $65^\circ$, and that the real cavity opening angle is not well constrained by their observations as it is degenerate with the inclination angle. The position angle of the Herbig–Haro objects observed in the near-IR varies between $258^\circ$ and $261^\circ$ (Poetzel et al. 1992), whilst Preibisch et al. (2003) adopted a value of $259^\circ$ for the outflow cavity symmetry axis even though the position angle of the loops ranges between $263^\circ$ and $265^\circ$ as determined from their K-band speckle interferometric observations.

The presence of a disc in AFGL 2591 is still uncertain. Through 2D axisymmetric radiative transfer modelling of the SED and 2MASS observations, Johnston et al. (2013) found that models with and without a disc can reproduce the observations that they modelled. However, HDQ, H$_2$$^{18}$O, and SO$_2$ interferometric observations at millimetre wavelengths point towards the presence of a sub-Keplerian disc-like structure and expanding material in the inner ~1000 au, and the continuum is extended nearly perpendicular to the jet/outflow direction (Wang, van der Tak & Hogerheijde 2012). The partially resolved source identified in the K-band speckle interferometric visibilities of Preibisch et al. (2003) seems to be tracing the inner rim of this disc, which has a size of ~120 au at 3.3 kpc as derived from their modelling.

The kinematics of the outflow and the inner region have been studied by several authors. The blue-shifted radial velocity of the jet is constrained between 200 and 500 km s$^{-1}$ as measured from the line wings of the Herbig–Haro objects (Poetzel et al. 1992) and

<footnote n="1">¹United Kingdom Infrared Telescope</footnote>

<footnote n="2">²James Clerk Maxwell Telescope</footnote>

<footnote n="3">³Atacama Pathfinder Experiment</footnote>

<footnote n="4">⁴The Herschel imaging survey of OB Young Stellar objects (HOBYS) is a Herschel key programme. See http://hobys-herschel.cea.fr</footnote>
<table><caption>**Table 1.** 70-µm angular sizes.</caption><thead><tr><th rowspan="2">Type</th><th colspan="2">2D Gaussian</th><th colspan="2">Horizontal slice</th><th colspan="2">Vertical slice</th></tr><tr><th>FWHM<br>(arcsec)</th><th>PA<br>(°)</th><th>FWHM<br>(arcsec)</th><th>FW1%<br>(arcsec)</th><th>FWHM<br>(arcsec)</th><th>FW1%<br>(arcsec)</th></tr></thead><tbody><tr><td>Observed</td><td>14.5 × 12.9</td><td>253.1 ± 0.7</td><td>12.2</td><td>58.0</td><td>11.3</td><td>47.7</td></tr><tr><td>Model</td><td>10.7 × 9.9</td><td>244 ± 0.8</td><td>8.2</td><td>48.9</td><td>8.0</td><td>42.6</td></tr><tr><td>PSF</td><td>5.9 × 6.6<sup>a</sup></td><td>56 ± 1</td><td>6.1</td><td>22.9</td><td>5.9</td><td>20.7</td></tr></tbody></table>

*Notes.* An error of 0.02 arcsec is estimated for the semi-axes of the 2D Gaussian fit. FW1% stands for full width at 1 per cent the peak intensity.

<sup>a</sup>Note this value is from the PSF binned to the observed pixel size and rotated to match the observations, hence it is different from the one presented in Table 5.

¹²CO (van der Tak et al. 1999). The C¹²O observations of Johnston et al. (2013) show evidence of Hubble-like expansion of the blue-shifted emission towards the outflow. Hubble-like expansion and rotation at disc scales (<5000 au) was also found from Plateau de Bure Interferometer (PdBI) observations of HDO and H₂¹²O by Wang et al. (2012). The rotation at these scales seems to be sub-Keplerian and anticlockwise looking down from the blue-shifted cavity. Wang et al. (2012) argue that these water isotopologues are closer to the surface of a disc-like structure, thus the expansion motion is a result of the interaction of the stellar radiation/wind, and infall occurs in the inner region of the disc where these molecules are depleted. Similarly, their SO₂ observations trace material affected more by the wind than the rotation.

The change in the abundance for different molecules in AFGL 2591 has been addressed by several authors. Benz et al. (2007) used the spherically symmetric abundance models of Stäuber et al. (2005) to fit their CS and SO interferometric (sub-arcsec resolution) observations, and concluded that X-ray emission from the protostar is needed to produce a better fit. Kaźmierczak-Barthel et al. (2015) studied the emission of 25 molecules observed with the Heterodyne Instrument for the Far-Infrared (HIFI) by Herschel (Pilbratt et al. 2010). They found that the emission from some molecules can be better described by abundances with a jump at 100 or 230 K, which were then fitted by a theoretical chemical model. Their theoretical model fitted relatively well the abundance jumps in molecules like H₂O and NH₃, and predicted a chemical age of 10–50 kyr. The first jump at 100 K is related to the evaporation of ices (van Dishoeck & Blake 1998) whilst the second jump at temperatures T > 230 K is important for N-bearing molecules due to less formation of N₂ (Boonman et al. 2001). Far Ultraviolet (UV) and X-ray radiation, which penetrates a thin layer in the cavity walls, can also enhance the abundance of several diatomic molecules observed in the far-IR and (sub)mm as shown by the chemical modelling of Bruderer et al. (2009), Bruderer et al. (2010). The physical–chemical modelling of 14 molecules observed at 1 mm with the Northern Extended Millimeter Array (NOEMA) by Gieser et al. (2019) shows that many species have abundances that change with radius (e.g. SO, SO₂), and derived a chemical age of ~20 kyr. They also obtained a flatter density distribution to explain the gas distribution (index of 1.0) than the one derived from dust emission at 1 mm (index of 1.7). These studies show that molecular line emission traces gas with different physical conditions, and imply the presence of different density structures from the changes in abundance as a function of temperature.

In this work, we have combined high-resolution observations at key frequencies to constrain the dust density, temperature, and gas velocity through radiative transfer modelling. The data are presented in Section 2. The dust continuum and molecular line modelling

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Observed (colour scale and grey contours) and model (green contours, see Section 3) 70-µm maps of AFGL 2591. Contour levels, in intensity normalized to peak value (87.3 Jy), are 5, 10, 20, ..., 1280 × σ<sub>rms</sub> with σ<sub>rms</sub> = 50 mJy. Bottom-hand and right-hand panels show a slice through the peak in the horizontal and vertical directions, respectively. The blue squares in the side panels correspond to the observations (errors are smaller than the point size), the green continuous line is the best-fitting model presented in Section 4 and the dotted cyan line is the PSF. The shaded region encloses the points higher than 5σ<sub>rms</sub>.</figcaption></figure>

procedures are presented in Section 3 and a detailed description of the procedures can be found in Appendices A–C (available online). The results of the modelling are presented in Section 4 and discussed in Section 5. Finally, our conclusions are presented in Section 6.

## 2 DATA

### 2.1 Herschel 70 µm

Resolved 70-µm data from Herschel observed as part of the HOBYS survey were used. The final images were obtained using procedures from the HOBYS team (Motte, private communication), and details of the reduction process can be found in Schneider et al. (2016). These data were taken with the Photodetector Array Camera and
<figure><img src="image_4.png" /><figcaption>**Figure 2.** NOEMA 1.3-mm observed (blue squares) and simulated (green circles, see Section 3) interferometry visibility profiles. (a) Compact (D) configuration visibility profile. (b) Extended (A+B) configurations visibility profile. The selection of the data to be fitted is explained in Appendix C (online).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** NOEMA 1.3-mm observed (colour scale and grey contours) and model (green contours, see Section 3) CLEAN images. The compact configuration map is shown in (a) and the levels are −2, 2, 3, 5, 10, 20 × σ<sub>rms</sub> with σ<sub>rms</sub> = 5.0 mJy beam<sup>−1</sup>. The white crosses mark the position of the continuum sources in Table 3. The orange star and circle show the position and 3.6-cm deconvolved size (~4.4 arcsec; Johnston et al. 2013) of VLA 1. The orange triangle shows the position of VLA 2. The extended configurations map is shown in (b) and the levels are −3, 3, 5, 10, 20, 40 × σ<sub>rms</sub> with σ<sub>rms</sub> = 0.8 mJy beam<sup>−1</sup>. For comparison, the VLA 3.6-cm observations of Johnston et al. (2013) are plotted in light blue contours at 3, 4, 5, 7, 10, 15, 20 × σ<sub>rms</sub> levels with σ<sub>rms</sub> = 30 µJy beam<sup>−1</sup>. The dark red square in (a) shows the area covered in (b). The observed beam is shown in the bottom left-hand corner.</figcaption></figure>

Spectrometer (PACS; Poglitsch et al. 2010) in the parallel mode of
the telescope at a scan speed of 20 arcsec s⁻¹. This speed produces
a more symmetrical PSF (6.51 × 5.44 arcsec²) with fewer artefacts
than data scanned at 60 arcsec s⁻¹, which has a 12.58 × 5.85 arcsec²
point spread function (PSF) (Lutz 2012; Olguin et al. 2015). A zoom
on AFGL 2591 is shown in Fig. 1 and provides the highest resolution
in the far-IR to date. Table 1 lists the observed size of the 70-µm
emission as obtained from a 2D Gaussian fit to the data, which shows
that the source is well resolved by Herschel and that the Gaussian
position angle is also close to the E–W outflow direction. The former
is also true for the horizontal and vertical slices shown in Fig. 1 side
panels. The emission peak is shifted ~3.2 arcsec to the west of the
radio source position (Trinidad et al. 2003). However, the Herschel
1σ astrometry error<fnref n="5" /> is ~2 arcsec, hence their positions agree within
2σ. In what follows, it will be assumed that the emission peak at
70-µm coincides with the radio one.

2.2 NOEMA 1.3-mm observations

AFGL 2591 was observed in the 1.3-mm spectral band with NOEMA
between 2014 August and 2016 February as part of the CORE
project (Beuther et al. 2018). The observations were undertaken in
three array configurations (A, B, and D) at 1.37 mm with the wide-
band correlator Widex and a narrow-band correlator (for a complete
overview of the spectral setup see Ahmadi et al. 2018). The data
were reduced by the CORE team and continuum visibilities were
extracted and then CLEANed (see Beuther et al. 2018). To facilitate
the comparison between models and observations, we separate the
analysis of the A+B visibilities, or extended configurations, and the
D visibilities, or compact configuration. For the purpose of this paper,
we use the continuum data from A+B and D configurations, and only
the compact configuration for the CH3CN molecular line data in order
to analyse the kinematics of the inner envelope (<10 000 au scales).
A study of the CH3CN chemistry using the combined A+B+D data
can be found in Gieser et al. (2019). The study of the kinematics
of the disc (if any) from the combined data set will be the focus of
future publications (e.g. Ahmadi et al. in preparation).

<footnote n="5">⁵http://herschel.esac.esa.int/Docs/PACS/pdf/pacs.om.pdf</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** NOEMA compact configuration zeroth (upper row) and first (lower row) moment maps of CH₃CN *J* = 12–11 for the *K* = 3 transition. The continuum is plotted in colour scale in the zeroth moment maps. The zeroth moment contour levels are 1, 3, 5, 10, 20, 40 × σ<sub>rms</sub> with σ<sub>rms</sub> = 0.45 Jy beam<sup>-1</sup> km s<sup>-1</sup>. Zero velocity corresponds to the systemic velocity in the LSR system (v<sub>LSR</sub> = -5.5 km s<sup>-1</sup>). The positions of the continuum sources are marked with a cross and the beam size is shown in the bottom left-hand corner.</figcaption></figure>

The beam size of the compact configuration continuum observations is 2.3 × 1.6 arcsec² PA = 887, whilst the beam size of the combined extended configurations is 0.48 × 0.38 arcsec² PA = 56°. Vector-averaged radial visibility profiles were extracted, where we used the standard error of the mean of the points in each uv-distance bin as measure of the errors. The profiles are shown in Fig. 2, whilst the CLEAN images are shown in Fig. 3. The NOEMA/CORE A+B observations recover more emission towards smaller baselines than previous PdBI 1.4-mm A+B observations by Wang et al. (2012) after taking into consideration the difference in wavelengths, whilst achieving similar angular resolution.

The methyl cyanide CH₃CN *J* = 12 – 11 line was covered by the narrow-band correlator with a spectral resolution of 0.312 MHz (0.43 km s⁻¹). A data cube covering the transition’s *K*-ladder between *K* = 0–5 was extracted. The beam size of these observations is 2.6 × 1.8 arcsec² PA = 87°. The average noise measured in an annulus of sky along all channels is σ<sub>rms</sub> = 30 mJy beam<sup>-1</sup> per channel, with variations between 10 and 50 mJy beam<sup>-1</sup> at different channels.

Zereth- and first-order moment maps for the transition with *K* = 3 are shown in Fig. 4. The first moment maps were

<table><caption>**Table 2.** Summary of NOEMA compact configuration continuum sources.</caption><thead><tr><th>NOEMA Source</th><th>VLA Source</th><th>RA (hh:mm:ss J2000)</th><th>Dec. (°:′:″ J2000)</th><th>Peak flux (mJy beam<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>3</td><td>20 29 24.864</td><td>40 11 19.51</td><td>230</td></tr><tr><td>2</td><td>1</td><td>20 29 24.627</td><td>40 11 15.22</td><td>62</td></tr><tr><td>3</td><td>–</td><td>20 29 24.900</td><td>40 11 15.60</td><td>61</td></tr></tbody></table>

calculated considering all data with fluxes larger than 5σ<sub>rms</sub> = 0.15 Jy beam<sup>-1</sup> per channel over a range of spectral channels covering the lines. The local standard of rest (LSR) velocity measured by Kaźmierczak-Barthel et al. (2014, -5.5 ± 0.5 km s<sup>-1</sup>) over 24 lines in Herschel/HIFI observations was used to shift the velocities, so zero velocity corresponds to the line systemic velocity.

We identified three continuum sources from the compact configuration data named NOEMA 1–3 and their positions and peak intensities are summarized in Table 2. The position of NOEMA 1 is consistent with the VLA 3 source, for consistency to previous works we will use the latter name. Fig. 3(a) reveals the presence of at least 2 sources south of VLA 3. Thus, the NOEMA images (continuum and zeroth moment compact configuration maps) were fitted with three 2D Gaussian components and the results are listed in Table 3. The observed size of VLA 3 decreases slightly with increasing *K* transition, i.e. decreases with increasing upper energy of the line, and the source is only partially resolved by the molecular line observations. The orientation of the source from the methyl cyanide emission seems to agree with the large-scale orientation of the near-IR images, but the beam is oriented in the same direction. In fact, the deconvolved PAs shown in Table 4 have large uncertainties, thus the orientation is strongly affected by the beam.

Similarly, a single 2D Gaussian was fitted to the extended configurations emission in Fig. 3(b). The results listed in Table 3 show that the position angle of the emission is closer to perpendicular to the radio jet in Fig. 5 (~100°, Johnston et al. 2013) rather than oriented in the outflow direction. Although the Gaussian fit is just an approximation, the orientation of the lower level contours in Fig. 3(b) shows a similar orientation.

To extract the kinematic information from the observations of VLA 3, we calculated position–velocity (pv) maps for CH₃CN lines with *K* = 2–5 and are shown in Fig. 6. The slices were obtained at 6° (along the putative disc major axis from Wang et al. 2012) and 276° (along the outflow direction; see online Appendix C for a discussion) with a slice width of 1.8 arcsec (three times the pixel size and close to the beamwidth).

The position of NOEMA 2 is close to VLA 1 and its total flux (~0.1 mJy) is consistent with the value predicted from the VLA 1 radio spectral index (α = 0, Johnston et al. 2013). In Fig. 5, it can be seen that the methyl cyanide emission of NOEMA 2 coincides with emission from the eastern hemisphere of VLA 1. This is consistent with VLA 1 being a slightly cometary H II region with the ionizing gas expanding towards the less dense gas in the west.

The source immediately south of VLA 3, NOEMA 3, has an intensity peak which is shifted ~4 arcsec towards the south-east from the continuum peak of VLA 3. This source is not associated with 3.6-cm radio continuum emission above a level of 3σ = 90 µJy beam<sup>-1</sup> (see Fig. 5) and has not been detected in the IR, thus its nature is as yet unknown. Fig. 4 also shows that the velocity structures of the three sources can be separated from each other.
<table><caption>**Table 3.** Measured properties from 2D Gaussian fit to the NOEMA 1.3-mm continuum and CH₃CN zeroth-order moment emission.</caption><thead><tr><th>K<sup>a</sup></th><th>Frequency (GHz)</th><th>E<sub>u</sub><sup>a</sup> (K)</th><th>RA (hh:mm:ss)</th><th>Dec. (°:′:′′)</th><th>Peak<sup>b</sup></th><th>T<sub>B</sub><sup>c</sup> (K)</th><th>Flux<sup>d</sup></th><th>Major (arcsec)</th><th>Minor (arcsec)</th><th>PA (°)</th></tr></thead><tbody><tr><td colspan="11" style="text-align:center; font-weight:bold;">NOEMA 1 (VLA 3)</td></tr><tr><td>Continuum (compact)</td><td></td><td></td><td>20 29 24.864</td><td>40 11 19.51</td><td>210 ± 6</td><td>-</td><td>436 ± 12</td><td>3.28 ± 0.04</td><td>2.37 ± 0.06</td><td>267 ± 2</td></tr><tr><td>Continuum (extended)</td><td></td><td></td><td>20 29 24.875</td><td>40 11 19.46</td><td>50 ± 3</td><td>-</td><td>180 ± 11</td><td>0.90 ± 0.05</td><td>0.72 ± 0.03</td><td>21 ± 9</td></tr><tr><td>0</td><td>220.74726</td><td>68.87</td><td>20 29 24.864</td><td>40 11 19.41</td><td>32.7 ± 0.6</td><td>40.9</td><td>39.3 ± 0.7</td><td>2.87 ± 0.03</td><td>2.04 ± 0.05</td><td>87.1 ± 0.9</td></tr><tr><td>1</td><td>220.74301</td><td>76.01</td><td>20 29 24.861</td><td>40 11 19.38</td><td>26.4 ± 0.5</td><td>37.9</td><td>31.6 ± 0.6</td><td>2.85 ± 0.03</td><td>2.04 ± 0.05</td><td>86.8 ± 1.0</td></tr><tr><td>2</td><td>220.73026</td><td>97.44</td><td>20 29 24.861</td><td>40 11 19.39</td><td>24.5 ± 0.6</td><td>34.5</td><td>28.9 ± 0.7</td><td>2.82 ± 0.04</td><td>2.03 ± 0.06</td><td>86.1 ± 1.3</td></tr><tr><td>3</td><td>220.70902</td><td>133.16</td><td>20 29 24.861</td><td>40 11 19.39</td><td>32.1 ± 0.6</td><td>53.9</td><td>37.8 ± 0.7</td><td>2.81 ± 0.04</td><td>2.03 ± 0.05</td><td>84.7 ± 1.1</td></tr><tr><td>4</td><td>220.67929</td><td>183.15</td><td>20 29 24.861</td><td>40 11 19.36</td><td>16.2 ± 0.4</td><td>22.3</td><td>18.3 ± 0.5</td><td>2.81 ± 0.05</td><td>1.96 ± 0.07</td><td>84.6 ± 1.3</td></tr><tr><td>5</td><td>220.64108</td><td>247.40</td><td>20 29 24.856</td><td>40 11 19.36</td><td>12.4 ± 0.3</td><td>15.3</td><td>13.7 ± 0.4</td><td>2.86 ± 0.05</td><td>1.93 ± 0.07</td><td>87.7 ± 1.3</td></tr><tr><td colspan="11" style="text-align:center; font-weight:bold;">NOEMA 2</td></tr><tr><td>Continuum (compact)</td><td></td><td></td><td>20 29 24.627</td><td>40 11 15.22</td><td>48 ± 6</td><td>-</td><td>76 ± 9</td><td>2.9 ± 0.3</td><td>2.0 ± 0.2</td><td>18 ± 7</td></tr><tr><td>0</td><td>220.74726</td><td>68.87</td><td>20 29 24.62</td><td>40 11 15.5</td><td>2.3 ± 0.3</td><td>2.7</td><td>3.8 ± 0.5</td><td>3.2 ± 0.3</td><td>2.6 ± 0.3</td><td>137 ± 14</td></tr><tr><td>1</td><td>220.74301</td><td>76.01</td><td>20 29 24.62</td><td>40 11 15.5</td><td>1.9 ± 0.3</td><td>3.3</td><td>3.3 ± 0.4</td><td>3.1 ± 0.3</td><td>2.7 ± 0.3</td><td>140 ± 15</td></tr><tr><td>2</td><td>220.73026</td><td>97.44</td><td>20 29 24.62</td><td>40 11 15.6</td><td>1.6 ± 0.3</td><td>2.3</td><td>3.4 ± 0.5</td><td>3.5 ± 0.4</td><td>3.0 ± 0.3</td><td>151 ± 17</td></tr><tr><td>3</td><td>220.70902</td><td>133.16</td><td>20 29 24.63</td><td>40 11 15.4</td><td>2.1 ± 0.3</td><td>2.9</td><td>4.1 ± 0.6</td><td>3.4 ± 0.4</td><td>2.8 ± 0.3</td><td>136 ± 15</td></tr><tr><td>4</td><td>220.67929</td><td>183.15</td><td>20 29 24.63</td><td>40 11 15.5</td><td>1.0 ± 0.2</td><td>1.3</td><td>2.1 ± 0.3</td><td>3.7 ± 0.4</td><td>2.8 ± 0.3</td><td>2.7 ± 9.5</td></tr><tr><td>5</td><td>220.64108</td><td>247.40</td><td>20 29 24.63</td><td>40 11 15.4</td><td>1.2 ± 0.1</td><td>1.1</td><td>2.0 ± 0.2</td><td>3.2 ± 0.2</td><td>2.6 ± 0.3</td><td>110 ± 9</td></tr><tr><td colspan="11" style="text-align:center; font-weight:bold;">NOEMA 3</td></tr><tr><td>Continuum (compact)</td><td></td><td></td><td>20 29 24.90</td><td>40 11 15.6</td><td>57 ± 6</td><td>-</td><td>250 ± 25</td><td>4.9 ± 0.2</td><td>3.3 ± 0.2</td><td>49 ± 5</td></tr><tr><td>0</td><td>220.74726</td><td>68.87</td><td>20 29 24.99</td><td>40 11 14.5</td><td>2.6 ± 0.3</td><td>2.8</td><td>3.2 ± 0.3</td><td>3.2 ± 0.2</td><td>1.9 ± 0.3</td><td>114 ± 4</td></tr><tr><td>1</td><td>220.74301</td><td>76.01</td><td>20 29 24.99</td><td>40 11 14.5</td><td>2.1 ± 0.2</td><td>3.1</td><td>2.8 ± 0.3</td><td>3.2 ± 0.2</td><td>2.0 ± 0.2</td><td>113 ± 4</td></tr><tr><td>2</td><td>220.73026</td><td>97.44</td><td>20 29 24.99</td><td>40 11 14.5</td><td>2.1 ± 0.3</td><td>2.9</td><td>2.3 ± 0.3</td><td>2.9 ± 0.3</td><td>1.8 ± 0.3</td><td>117 ± 6</td></tr><tr><td>3</td><td>220.70902</td><td>133.16</td><td>20 29 24.98</td><td>40 11 14.5</td><td>2.5 ± 0.3</td><td>3.4</td><td>3.2 ± 0.4</td><td>3.3 ± 0.2</td><td>1.9 ± 0.3</td><td>110 ± 4</td></tr><tr><td>4</td><td>220.67929</td><td>183.15</td><td>20 29 24.98</td><td>40 11 14.4</td><td>1.5 ± 0.3</td><td>1.4</td><td>1.7 ± 0.3</td><td>3.0 ± 0.3</td><td>1.8 ± 0.4</td><td>111 ± 6</td></tr><tr><td>5</td><td>220.64108</td><td>247.40</td><td>20 29 24.97</td><td>40 22 14.6</td><td>1.0 ± 0.2</td><td>1.3</td><td>1.2 ± 0.2</td><td>3.2 ± 0.3</td><td>1.8 ± 0.4</td><td>110 ± 6</td></tr></tbody></table>

The 2D Gaussian components were fitted simultaneously to the continuum images, whilst for the methyl cyanide zeroth moment maps each source was fitted separately within a box isolating each source.

<table><caption>**Table 4.** Deconvolved 2D Gaussian fit sizes for NOEMA 1 (VLA 3) from NOEMA compact configuration CH₃CN emission at different *K* transition.</caption><thead><tr><th>K</th><th>Major (arcsec)</th><th>Minor (arcsec)</th><th>PA (°)</th></tr></thead><tbody><tr><td>0</td><td>1.14 ± 0.08</td><td>0.86 ± 0.12</td><td>87 ± 175</td></tr><tr><td>1</td><td>1.09 ± 0.09</td><td>0.88 ± 0.13</td><td>83 ± 173</td></tr><tr><td>2</td><td>1.00 ± 0.13</td><td>0.84 ± 0.19</td><td>73 ± 79</td></tr><tr><td>3</td><td>1.01 ± 0.13</td><td>0.80 ± 0.18</td><td>58 ± 41</td></tr><tr><td>4</td><td>1.00 ± 0.17</td><td>0.61 ± 0.37</td><td>69 ± 31</td></tr><tr><td>5</td><td>1.12 ± 0.14</td><td>0.55 ± 0.43</td><td>90 ± 176</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** VLA 3.6-cm continuum emission (green contours) from Johnston et al. (2013) superimposed to the NOEMA compact configuration CH₃CN *J* = 12 − 11 *K* = 3 zeroth moment map (colour scale). Sources are labelled and the position of the continuum source NOEMA 3 is marked with a white star. Following Johnston et al. (2013), contours are 3, 4, 5, 7, 10, 15, 20, 30–100 × σ<sub>rms</sub> (σ<sub>rms</sub> = 30 µJy beam⁻¹).</figcaption></figure>

## 2.3 Other continuum multiwavelength data

Multiwavelength observations were used to constrain the dust density and temperature distributions. Since these images were previously processed, only the sky level was subtracted when necessary. Table 5 presents a summary of the observations, which were used to extract spatial information for the fitting and are described below.

### 2.3.1 Near-IR imaging

Near-IR images from the United Kingdom Infrared Telescope (UKIRT) Wide Field Camera (WFCAM, Casali et al. 2007) at 1.2, 1.6, and 2.1 µm, *J*, *H* and *K* bands, respectively, were obtained from the UKIRT Infrared Deep Sky Survey (UKIDSS; Lucas et al. 2008) and are presented in Fig. 7. These data have higher resolution

(~1 arcsec) than the 2MASS data used by Johnston et al. (2013), which had a resolution of ~2.5 arcsec. The images were flux-calibrated and error maps were obtained by using Poisson statistics on the un-calibrated data converted to counts. A Moffat function was

<footnote>ᵃFrom Leiden Atomic and Molecular Database (LAMDA).</footnote>

<footnote>ᵇThe units of peak fluxes are mJy beam⁻¹ for the continuum observations and Jy beam⁻¹ km s⁻¹ for the molecular lines.</footnote>

<footnote>ᶜBrightness temperature at the line peak at the position of the zeroth moment peak.</footnote>

<footnote>ᵈThe units of flux are mJy for the continuum observations and Jy km s⁻¹ for the molecular lines.</footnote>
<table><caption>**Table 5.** Summary of continuum observations used to extract spatial information.</caption><thead><tr><th>λ<br>(μm)</th><th>Instrument</th><th>Resolution</th><th>Type<sup>a</sup></th></tr></thead><tbody><tr><td>1.2</td><td>UKIRT/WFCAM</td><td>0.9 arcsec</td><td>I</td></tr><tr><td>1.6</td><td>UKIRT/WFCAM</td><td>0.9 arcsec</td><td>I</td></tr><tr><td>2.1</td><td>UKIRT/WFCAM</td><td>0.9 arcsec</td><td>I</td></tr><tr><td>2.1</td><td>SAO/6 m</td><td>0.17 arcsec</td><td>V</td></tr><tr><td>70</td><td>Herschel/PACS</td><td>6.51 × 5.442 arcsec<sup>2</sup></td><td>I</td></tr><tr><td>450</td><td>JCMT/SCUBA</td><td>9 arcsec</td><td>P</td></tr><tr><td>850</td><td>JCMT/SCUBA</td><td>14 arcsec</td><td>P</td></tr><tr><td>1300</td><td>NOEMA (D)</td><td>2.3 × 1.6 arcsec<sup>2</sup> PA 887</td><td>V</td></tr><tr><td>1300</td><td>NOEMA (A+B)</td><td>0.48 × 0.38 arcsec<sup>2</sup> PA 56°</td><td>V</td></tr></tbody></table>

fitted to several saturated and unsaturated nearby point-like sources to produce a PSF image. This function has been proven reliable in reproducing the PSF wings (e.g. McDonald et al. 2011). The core of saturated sources was masked in order to fit better the PSF wings. The full width at half-maximum (FWHM) of the fitted Moffat function are 0.9 ± 0.1, 0.9 ± 0.2, and 0.9 ± 0.2 arcsec and the atmospheric scattering coefficients are 3.0 ± 0.4, 2.7 ± 0.4, and 2.7 ± 0.7 for the *J*, *H*, and *K* bands, respectively. These widths are consistent with the observed seeing FWHM of 0.8 arcsec, as recorded in the header of the images.

In order to study the regions close to the MYSO, which are saturated in the UKIRT images, we used Special Astrophysical Observatory (SAO) 6-m telescope *K*-band speckle interferometric visibilities from Preibisch et al. (2003). The visibilities were averaged over annuli with constant width and the standard deviation of the data in each annulus was used as a measure of the errors. Fig. 8 shows the visibility radial profile. It is worth noticing that the visibilities are reasonably symmetric, which results in relatively small errors (see Preibisch et al. 2003, their fig. 2).

### 2.3.2 (Sub)mm

In the submm, images were obtained from the JCMT Legacy Catalogue (Di Francesco et al. 2008), which were taken with SCUBA<fnref n="6" /> at 450 and 850 µm. Following Di Francesco et al. (2008), we prepared a PSF composed of two Gaussians, which reproduce the main and secondary lobes. The main lobes have an FWHM of 9 and 14 arcsec at 450 and 850 µm, respectively, whilst 2D Gaussians fitted to observed emission towards AFGL 2591 have an FWHM of 37.1 × 28.5 arcsec² PA = 247° ± 1° and 47.6 × 38.7 arcsec² PA = 233° ± 1°, respectively. Hence, the source is elongated with position angles of the 2D Gaussian major axis consistent with the outflow cavity direction, as in the 70-µm data. However, since the real PSF is known to be more complicated and changes between nights (e.g. Hatchell et al. 2000), we fit azimuthally averaged radial intensity profiles (hereafter radial intensity profiles). The radial intensity profiles shown in Fig. 9 were calculated as the average flux of concentric annuli regions of constant width and their errors are the standard deviation of the enclosed fluxes.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Observed pv maps (grey scale and blue contours) and best-fitting model pv maps (green contours, see Section 3) of the CH₃CN *J* = 12–11, *K* = 2–5 line emission for a position angle of 276° (left-hand column), i.e. along the rotation axis, and 6° (right-hand column). Contours levels are 10–80 per cent in 20 per cent steps of 5, 8, 3, and 2 Jy beam⁻¹ for each *K* line, respectively.</figcaption></figure>

<footnote><sup>a</sup>Type of data product used to compare observations and models: I for images, V for visibility profiles, and P for azimuthally averaged radial profiles.</footnote>

<footnote n="6">⁶Submillimetre Common-User Bolometer Array</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Observed (colour map and grey contours) and model (green contours, see Section 3) near-IR maps of AFGL 2591. The J-band data and model are shown in (a) and contour levels are 3, 5, 10, 20, 40 × σ<sub>rms</sub> with σ<sub>rms</sub> = 0.3 MJy sr<sup>-1</sup>. (b) shows the H-band data and model in contours at levels 5, 10, 20–160 × σ<sub>rms</sub> with σ<sub>rms</sub> = 0.4 MJy sr<sup>-1</sup>. (c) shows the K-band data and model in contours at levels 5, 10, 20–2560 × σ<sub>rms</sub> with σ<sub>rms</sub> = 0.6 MJy sr<sup>-1</sup>. The colour maps are in units of MJy sr<sup>-1</sup>.</figcaption></figure>

2.3.3 *Spectral energy distribution*

The SED consists of data points from Johnston et al. (2013) for
$\lambda < 60$ µm and the data points in Table 6 for longer wavelengths.
The whole SED is shown in Fig. 10. The fluxes for the Herschel
bands, observed by PACS, and the Spectral and Photometric Imaging

<figure><img src="image_10.png" /><figcaption>**Figure 8.** K-band speckle interferometry visibilities radial profile. Observed points are represented by blue squares and the best-fitting model by a continuous green line. The model clearly underestimates the observations (see Section 4.1).</figcaption></figure>

Receiver (SPIRE; Griffin et al. 2010), were obtained from the sum of all pixels with intensity higher than three times the image noise. This filters out the filamentary emission observed in SPIRE data whilst the fluxes are within the calibration uncertainties of the ones calculated with an aperture covering the whole region. The flux calibration uncertainties for PACS and SPIRE in the literature range from ~5 to 15 per cent (e.g. Sadavoy et al. 2013; Roy et al. 2014, and Herschel online documentation^7). We adopted an uncertainty of 10 per cent for both.

3 DATA MODELLING

**3.1 Multiwavelength data modelling**

We performed a radiative transfer modelling of the multiwave-
length dust continuum and CH3CN observations to obtain density,
temperature, and velocity distributions. The fit was performed in
three steps and is described in Fig. 11(a). Each step involved the
production of synthetic observations as described in Fig. 11(b) to
obtain the model data products (images, intensity radial profiles,
etc.) which were then compared with the observations. To produce
synthetic data, we used the 3D radiative transfer codes HYPERION
(Robitaille 2011) for dust continuum and MOLLIE (Keto & Rybicki
2010) for line emission. The modelling steps included one visual
inspection for the continuum fitting and model grid fitting for some
continuum observations and the CH3CN observation. Details and
the particularities of the modelling and synthetic observations are
described in Appendices A–C (online).

## 3.2 1D local thermodynamic equilibrium (LTE) and non-LTE CH3CN modelling

We fitted the methyl cyanide 12–11 *K*-ladder with the program
CASSIS<fnref n="8" /> to study some physical properties of the gas (e.g. tempera-

<footnote n="7">⁷SPIRE: http://herschel.esac.esa.int/twiki/bin/view/Public/SpireCalibration
Web</footnote>

<footnote>PACS: http://herschel.esac.esa.int/twiki/bin/view/Public/PacsCalibrationW
eb</footnote>

<footnote n="8">⁸Based on analysis carried out with the CASSIS software and the Cologne Database for Molecular Spectroscopy (CDMS; https://cdms.astro.uni-koeln.de) molecular database. CASSIS has been developed by IRAP-UPS/CNRS (http://cassis.irap.omp.eu).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Normalized radial intensity profiles at 450 and 850 µm. Observed profiles are represented with blue squares and the best-fitting model by a green continuous line (see Section 3). The horizontal dotted line represents the $3\sigma$ level above zero, where $\sigma$ is the noise level. Observed points above this level were compared with the model ones. The cyan dotted line shows the PSF profile.</figcaption></figure>

<table><caption>**Table 6.** SED data points for $\lambda \ge 60$ µm.</caption><thead><tr><th>Wavelength (µm)</th><th>Flux (Jy)</th><th>Instrument</th><th>Reference</th></tr></thead><tbody><tr><td>70</td><td>5600 ± 560</td><td>Herschel/PACS</td><td>(1)</td></tr><tr><td>160</td><td>4055 ± 406</td><td>Herschel/PACS</td><td>(1)</td></tr><tr><td>250</td><td>1250 ± 125</td><td>Herschel/SPIRE</td><td>(1)</td></tr><tr><td>350</td><td>426 ± 43</td><td>Herschel/SPIRE</td><td>(1)</td></tr><tr><td>450</td><td>222 ± 110</td><td>JCMT/SCUBA</td><td>(2)</td></tr><tr><td>500</td><td>123 ± 12</td><td>Herschel/SPIRE</td><td>(1)</td></tr><tr><td>850</td><td>26.8 ± 5.2</td><td>JCMT/SCUBA</td><td>(2)</td></tr><tr><td>1200</td><td>8.1 ± 0.8</td><td>IRAM<br>30m/MAMBO</td><td>(3)</td></tr></tbody></table>

(1) – This work. Data from PACS and SPIRE were obtained from the HOBYS observations, with SPIRE data from high gain dedicated observations.

(2) – Di Francesco et al. (2008).

(3) – Roy et al. (2011) based on Motte et al. (2007).

ture). This program uses a Markov chain Monte Carlo, which samples the parameter space in order to minimize the $\chi^2$ between the observed and synthetic spectra. We used one isothermal density component model to fit the observations and their parameter ranges are listed in Table 7. CASSIS calculates spectra under the LTE approximation, or the non-LTE approximation by using the radiative transfer code RADEX (van der Tak et al. 2007). Two-density component models were experimented upon, but similar results were obtained.

In order to limit the region to analyse, only spectra with at least one line with a flux higher than $0.25$ Jy beam$^{-1}$ ($>5\sigma$) were considered in the fit. Additionally, CASSIS requires an estimate of each line rms and the flux calibration error. The first was estimated from boxes around the lines and vary between 0.04 and 0.06 K, where the higher the K line the higher the error. The relative flux calibration error was fixed to 20 per cent (Beuther et al. 2018).

The size parameter in Table 7 is the source size and is used to calculate the beam dilution factor, which is defined in CASSIS as

$$ \Omega = \frac{\theta^2}{\theta^2 + \theta_B^2}, \qquad (1) $$

where $\theta$ is the source size and $\theta_B$ is the beam size derived from the configuration files specific for each telescope, which include its diameter/largest baseline. In this case, the baseline used was

<figure><img src="image_12.png" /><figcaption>**Figure 10.** SED of AFGL 2591 and the best-fitting model. The emission from each of the different components included in the model are also shown (see legend).</figcaption></figure>

<table><caption>**Table 7.** Parameter constraints of the modelling with CASSIS.</caption><thead><tr><th>Parameter</th><th>Limits</th></tr></thead><tbody><tr><td>Column density ($N_{CH_3CN}$)</td><td>10<sup>13</sup>–10<sup>18</sup> cm<sup>-2</sup></td></tr><tr><td>H<sub>2</sub> density ($n_{H_2}$)<sup>a</sup></td><td>10<sup>3</sup>–10<sup>9</sup> cm<sup>-3</sup></td></tr><tr><td>Excitation temperature ($T_{ex}$)</td><td>10–500 K</td></tr><tr><td>Size</td><td>0–1.2 arcsec</td></tr><tr><td>Line width ($\Delta v$)</td><td>1–7 km s<sup>-1</sup></td></tr><tr><td>Offset to LSR velocity ($\Delta v_{LSR}$)<sup>b</sup></td><td>-2–2 km s<sup>-1</sup></td></tr></tbody></table>

<sup>a</sup>Only for the non-LTE case.

<sup>b</sup>$v_{LSR} = -5.5$ km s$^{-1}$.

calculated to match the resolution of the observations. The brightness temperature is then calculated by CASSIS as

$$ T_b = \Omega J_v(T_{ex})(1 - e^{-\tau}) + f(\text{CMB, dust}), \qquad (2) $$
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Radiative transfer modelling procedure. (a) Overall procedure. The blue dashed boxes show the two steps followed to fit the dust continuum data, whilst the grey dashed box shows the procedure used to fit the CH₃CN molecular line emission observations. Purple circles involve performing radiative transfer steps. (b) Image processing procedure performed at the end of each radiative transfer computation. Green boxes represent final products of each step. Detail descriptions of the procedures and additional tables are available in the online Appendices.</figcaption></figure>
<table><caption>**Table 8.** Parameters of the continuum and CH$_3$CN line best-fitting model, which consists of an Ulrich envelope and a flared disc.</caption><thead><tr><th>Parameter</th><th>Visual inspection</th><th>Grid</th><th>Constrained by</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center;">Continuum</td></tr><tr><td>Distance (kpc)</td><td></td><td>3.3 ± 0.1</td><td>Fixed</td></tr><tr><td>Position angle (PA)</td><td></td><td>259°</td><td>Fixed</td></tr><tr><td>Inclination angle (i)</td><td>40°</td><td>25° ± 5°</td><td>SED + near-IR</td></tr><tr><td>Stellar temperature (T<sub>*</sub>/K)</td><td>20000</td><td>18000</td><td>SED</td></tr><tr><td>Stellar luminosity (L<sub>*</sub>/10<sup>5</sup> L<sub>☉</sub>)</td><td>1.6</td><td>1.6 ± 0.2</td><td>SED</td></tr><tr><td>Inner radius (R<sub>in</sub>/au)</td><td>36</td><td>35</td><td>SED + Speckle</td></tr><tr><td>Outer radius (R<sub>out</sub>/au)</td><td>195000</td><td>195000</td><td>Submm</td></tr><tr><td>Envelope infall rate (M&#775;<sub>env</sub>/10<sup>-3</sup> M<sub>☉</sub> yr<sup>-1</sup>)</td><td>1.1</td><td>2.0 ± 0.2</td><td>SED + (sub)mm</td></tr><tr><td>Envelope dust</td><td></td><td>OHM92</td><td>NOEMA ext.</td></tr><tr><td>Centrifugal radius (R<sub>c</sub>/au)</td><td>200</td><td>2200 ± 200</td><td>NOEMA ext.</td></tr><tr><td>Disc mass (M&#775;<sub>d</sub>/M<sub>☉</sub>)</td><td>1.0</td><td>6 ± 1</td><td>NOEMA ext.</td></tr><tr><td>Disc scale height at r = 100 au (h<sub>0</sub>/au)</td><td></td><td>0.6</td><td>NOEMA ext.</td></tr><tr><td>Disc flaring exponent (α)</td><td>1.88</td><td>2.25</td><td>NOEMA ext.</td></tr><tr><td>Disc vertical density exponent (β)</td><td>1.13</td><td>1.25</td><td>NOEMA ext.</td></tr><tr><td>Disc dust</td><td></td><td>OHM92</td><td>Near-IR + SED</td></tr><tr><td>Cavity opening angle<sup>a</sup> (2θ<sub>cav</sub>)</td><td>60°</td><td>60° ± 5°</td><td>Near-IR + SED</td></tr><tr><td>Cavity shape exponent (b<sub>cav</sub>)</td><td></td><td>2.0</td><td>Near-IR + SED</td></tr><tr><td>Cavity density exponent (p<sub>cav</sub>)</td><td></td><td>2.0</td><td>Near-IR + SED</td></tr><tr><td>Cavity reference density<sup>a</sup> (ρ<sub>cav</sub>/10<sup>-21</sup> g cm<sup>-3</sup>)</td><td></td><td>1.4</td><td>Near-IR + SED</td></tr><tr><td>Cavity dust</td><td></td><td>KMH</td><td>Near-IR + SED</td></tr><tr><td colspan="4" style="text-align: center;">Extra parameters for line emission</td></tr><tr><td>Position angle (PA)</td><td>—</td><td>276°</td><td>Fixed</td></tr><tr><td>Stellar mass (M<sub>*</sub>/M<sub>☉</sub>)</td><td>—</td><td>35</td><td>Fixed</td></tr><tr><td>Abundance for T > 100 K (X<sub>1</sub>)</td><td>—</td><td>2 × 10<sup>-8</sup></td><td>Fixed</td></tr><tr><td>Turbulent width (Δv<sub>nth</sub>/km s<sup>-1</sup>)</td><td>—</td><td>2</td><td>Fixed</td></tr><tr><td>Outflow velocity<sup>a</sup> (v<sub>0</sub>/km s<sup>-1</sup>)</td><td>—</td><td>0, 1<sup>b</sup></td><td>Fixed</td></tr></tbody></table>

Parameters in bold under the last column are the 10 free parameters that can be constrained by the observations. The inner radius depends on the stellar luminosity.

<sup>a</sup>Defined at a height z = 10<sup>4</sup> au.

<sup>b</sup>2 models had the same best overall ranking value.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Unprocessed model image of the best-fitting model at 70 µm. Contour levels are 0.1, 0.3, 1, and 10 per cent of the peak intensity. Note the elongation along the cavity axis oriented at PA = 259° as the distance to the source increases.</figcaption></figure>

where $J_v(T_{ex}) = h\nu k^{-1}(e^{hv/kT_{ex}} - 1)^{-1}$, $T_{ex}$ is the excitation temperature, $\tau$ is the line optical depth and f(CMB, dust) accounts for the contribution from the cosmic microwave background and dust emission.

# 4 RESULTS

## 4.1 Continuum modelling

The parameters of the overall best-fitting model, i.e. after the continuum grid fitting step, are listed in Table 8. The extended emission observed at 70 µm is relatively well fitted by the best-fitting model up to ~40 arcsec from the centre as shown in Fig. 1. In general, the model matches the elongation in the data. At 1 per cent of the peak intensity, where the model has a better fit, the horizontal slice, which cuts through the cavity, is more elongated than the vertical one in the model and the observations as listed in Table 1. This cannot be explained by any asymmetry in the 70-µm PSF, which is relatively symmetric as shown in the same table. Fig. 12 shows the 70-µm image before convolution with the PSF, where the contribution of the cavity to the elongation can be observed. The emission is underestimated close to the source and along the cavity axis (PA = 259°, see Fig. 1). As a result, the angular sizes of the model at the FWHM from slices and 2D Gaussian fit listed in Table 1 are similar between each direction.

The parameters of the best-fitting model during the visual inspection step are listed in Table 8 and the ranking and reduced $\chi^2$ values for each observation are listed in Table 9. The best-fitting model has a type A density distribution (Ulrich 1976 envelope and a flared disc). The envelope radius was well-fitted during this step and has a value roughly a factor of 1.4 smaller than the initial visual inspection model and has a similar size to the model with disc of Johnston et al. (2013). The disc scale height is smaller by roughly one order
<table><caption>**Table 9.** Overall ranking and reduced $\chi^2$ values for the continuum best-fitting models of each type.</caption><thead><tr><th rowspan="2">Type</th><th rowspan="2">Rank<sup>a</sup></th><th rowspan="2">70 µm</th><th rowspan="2">SED</th><th rowspan="2">J</th><th rowspan="2">H</th><th colspan="2">K</th><th rowspan="2">450 µm</th><th rowspan="2">850 µm</th><th colspan="2">1.3 mm</th><th rowspan="2">χ<sup>2</sup></th><th rowspan="2">⟨χ<sup>2</sup>⟩</th></tr><tr><th>UKIRT</th><th>Speckle</th><th>Extended</th><th>Compact</th></tr></thead><tbody><tr><td colspan="14" style="text-align:center; font-weight:bold;">Overall best-fitting model</td></tr><tr><td>A</td><td>–</td><td>211</td><td>7.0</td><td>26</td><td>64</td><td>316</td><td>32</td><td>13</td><td>0.40</td><td>408</td><td>557</td><td>–</td><td>–</td></tr><tr><td colspan="14" style="text-align:center; font-weight:bold;">Best-fitting model by density type</td></tr><tr><td>A</td><td>662</td><td>114</td><td>9.7</td><td>23</td><td>79</td><td>311</td><td>49</td><td>20</td><td>0.56</td><td>993</td><td>852</td><td>137</td><td>164</td></tr><tr><td>B</td><td>854</td><td>160</td><td>7.3</td><td>25</td><td>81</td><td>305</td><td>35</td><td>17</td><td>0.51</td><td>1040</td><td>1142</td><td>137</td><td>194</td></tr><tr><td>C</td><td>1019</td><td>305</td><td>44.1</td><td>21</td><td>107</td><td>523</td><td>469</td><td>3</td><td>0.10</td><td>557</td><td>300</td><td>217</td><td>175</td></tr><tr><td colspan="14" style="text-align:center; font-weight:bold;">Best-fitting model by χ<sup>2</sup></td></tr><tr><td>A</td><td>1134</td><td>100</td><td>8.9</td><td>26</td><td>91</td><td>223</td><td>44</td><td>45</td><td>3.60</td><td>926</td><td>1570</td><td>113</td><td>160</td></tr><tr><td colspan="14" style="text-align:center; font-weight:bold;">Best-fitting model by ⟨χ<sup>2</sup>⟩</td></tr><tr><td>C</td><td>1278</td><td>274</td><td>54.2</td><td>23</td><td>155</td><td>624</td><td>25</td><td>0.7</td><td>1.16</td><td>156</td><td>195</td><td>267</td><td>125</td></tr></tbody></table>

of magnitude than those in Johnston et al. (2013), but this parameter cannot be constrained by our observations, thus it was less explored. The shape of the cavity changed from the initial exponent of 1.5 to 2.0, and its reference density decreased an ~15 per cent from its initial value and its density exponent did not change.

Similarly, in Table 8 and 9, we list the results for the overall best-fitting model after the grid-fitting step. The disc mass and radius increased with respect to the visual inspection best-fitting model. The former is roughly a third of the disc mass in the control model of Johnston et al. (2013), but the disc radii are similar. This difference is explained by the difference in the disc scale height. The cavity opening angle is equal to the one from the visual inspection best-fitting model, and the inclination angles differ by 15°. The luminosities are the same between the best-fitting models at each step, whilst the envelope infall rate is ~1.8 times higher in the overall best-fitting model.

The model SED agrees reasonably well with the observations as shown in Fig. 10. The total mass of the model, $M = 2.1 \times 10^3 M_\odot$, is enough to fit the 10-µm silicate feature relatively well at an inclination of 25° and closely fit the submm points. This is 900 $M_\odot$ more massive than the model with disc of Johnston et al. (2013), but 600 $M_\odot$ less massive than their envelope only model. Models with WD01 dust can also fit in the submm with lower masses, however, they underestimate the fluxes between 3 and 8 µm by a factor of ~10.

The near-IR maps in Fig. 7(a)–(c) show that the near-IR band images are well reproduced. However, Fig. 8 shows that the small scales closer to the star position probed by the speckle observation are not well reproduced. In general, most models with a cavity density of the same order as the best-fitting model do not fit this observation.

At 450 µm, the best-fitting model overestimates the emission at larger scales as shown in Fig. 9. Type C models with WD01 dust are able to reproduce this profile but not the near-IR observations. At 850 µm, most of the models reproduce the radial profile relatively well. These two fits are unaffected by the inclination angle. Palau et al. (2014) successfully fitted these profiles with a spherically symmetric model and obtained a density power-law index of 1.8, which is slightly lower than our type C models but steeper than the Ulrich envelope. However, an analysis of the observed images shows that the observed source is more elongated along the cavity at both wavelengths, but this is not seen in our models.

The 1.3-mm extended configuration visibilities are relatively well fitted as shown in Fig. 2(b). Models with WD01 dust or larger and more massive discs are able to match the data for baselines smaller than ~150 kλ, but these baselines were not included in the fitting process. For longer baselines, the intensity is determined mainly by

the mass of the disc, $6 M_\odot$ in the best-fitting model. The best-fitting model underestimates the observed 1.3-mm compact configuration visibilities (Fig. 2a) at large spatial scales (smaller uv distances). This is observed in most models. However, those models with WD01 dust give an improvement of the fit, but they do not fit the short wavelengths (λ < 12 µm) of the SED.

At ~10⁴ au scales, as mapped by the 1.3-mm compact configuration observations in Fig. 3(a), the model FWHM is $3.2 \times 3.0 \text{arcsec}^2$ PA = 325° ± 2° as obtained from a 2D Gaussian fit. This is slightly more extended than the observed size in Table 3 and has a different orientation. The presence of at least one source to the south of VLA 3, as shown in Fig. 3(a), may also contaminate the submm observations making them look more elongated than they really are. A more complete modelling of VLA 3 including the contribution of NOEMA 1 is beyond the scope of this paper.

Fig. 3(b) shows the degree of agreement between the observations and model at ~100 au scales from the 1.3-mm extended configuration observations. It shows that the model emission is more compact than the observed one, with an FWHM of $0.7 \times 0.6 \text{arcsec}^2$ PA = 247° ± 10° from a 2D Gaussian fit (cf. Table 3). The PA of the model is closer to the orientation of the outflow cavity than to the orientation of the disc, which can be explained by the inclination angle. The extended emission in the observations also show features, which cannot be explained by the model due to its symmetry.

### 4.1.1 Other density distributions

Table 9 lists the $\chi^2$ values for the best-fitting models for the different types of model densities (A, B, and C) for each observation in the visual inspection fitting. The values of the parameters of the best-fitting models are listed in Table D1 (online). The best-fitting type B model has the same stellar properties as the best-fitting type A model (from grid fitting). The envelope of the best-fitting type B model has the same radius and has an infall rate $0.5 \times 10^{-3} M_\odot \text{yr}^{-1}$ lower than the best-fitting type A model. The inclination angle is 10° lower than the best-fitting type A model, which may explain the better fit to the SED than the best-fitting type A from visual inspection. The distribution and shape of the cavity is the same as the best-fitting type A model. This model has an accretion luminosity of ~20 $L_\odot$.

The best-fitting type C model has the same inclination angle as the best-fitting type A model. It fits relatively well the submm and mm observations, but is worse at fitting the near-IR observations and the SED because its envelope has a WD01 dust. The envelope has a mass of $1.0 \times 10^3 M_\odot$ and the density in the cavity is lower ($\rho_{0,cav} = 3.5 \times 10^{-24} \text{ g cm}^{-3}$) but decreases at the same rate as the

<footnote><sup>a</sup>Average weighted ranking as described in Appendix A1.2 (online).</footnote>
best-fitting type A model. Its disc is more compact than the best-fitting type A model one from visual inspection, and slightly more massive (1.2 M☉). In general, this type of model does not fit the 70-μm image well.

## 4.2 CH$_3$CN line modelling

### 4.2.1 1D CASSIS modelling

We obtained similar fits to the data and results under the LTE and non-LTE approximations. The results for the LTE modelling are shown in Fig. 13. The peak column density for VLA 3 is $4.2 \times 10^{16} \text{ cm}^{-2}$ in the LTE model, and is located at the position of the source. The average column density inside a region with radius 2 arcsec is $1.5 \times 10^{16} \text{ cm}^{-2}$. An approximate value of the H$_2$ column density can be obtained from the continuum map in the optically thin limit (e.g. Morales et al. 2009):

$$N_{\text{H}_2} = \frac{S_v R_{\text{gd}}}{\Omega_B B_v(T_d) \kappa_v \mu m_{\text{H}}} \quad (3)$$

with $S_v$ the observed dust continuum flux and $\Omega_B$ the solid angle of the beam. If equation (3) is solved with a dust temperature $T_d = 190$ K, a dust opacity $\kappa_{1.3\text{mm}} = 1 \text{ cm}^2 \text{ g}^{-1}$, a mean molecular weight $\mu = 2.3$ and a gas-to-dust ratio $R_{\text{gd}} = 100$, a peak H$_2$ column density value of $2.0 \times 10^{23} \text{ cm}^{-2}$, and an average of $9.9 \times 10^{22} \text{ cm}^{-2}$ within the same 2-arcsec region are obtained. Hence, the methyl cyanide abundance is $2.1 \times 10^{-7}$ for the peak value and $1.5 \times 10^{-7}$ for the average.

The excitation temperature for the LTE model in VLA 3 peaks at the source position and has a value of 271 K. In a region of radius 2 arcsec, the minimum temperature is 120 K and the average is 190 K.

The FWHM of the line gets wider along the blue-shifted outflow direction. This may be due to turbulent motion of the gas in the outflow. The redshifted outflow direction does not show the same trend, probably because the emission is optically thick towards that side. However, the fit is worse towards the red-shifted outflow direction, as shown by the reduced $\chi^2$ maps, hence this result should be taken with caution. The line velocity maps of VLA 3 are consistent with the first moment maps in Fig. 4.

The sizes obtained from the modelling are smaller than the beam size of the observations and the deconvolved angular size derived from continuum observations (~2.0 arcsec from Table 3). The less beam diluted areas are located towards the denser regions. As the flux is proportional to the beam dilution factor $\Omega$ and proportional to the column density in the optically thin regime, the size and the column density parameters may be degenerate. However, if the core was spherical, the depth along the line of sight would increase towards the centre of the core.

The regions NOEMA 2 and 3 show a different temperature and velocity structure. As expected for an H II region (VLA 1), NOEMA 2 has wider lines due to turbulent gas expansion and is hotter than the other two regions. NOEMA 2 also has a lower column density in comparison to the other two regions, which is consistent with it being optically thin at mm/radio wavelengths as found by e.g. Trinidad et al. (2003) at 3.6 cm. NOEMA 3 peaks in temperature and column density at the same position. Its peak temperature is 240 K and its column density is $1.6 \times 10^{16} \text{ cm}^{-2}$, hence is slightly colder and less dense than VLA 3. There appears to be a similar velocity gradient in NOEMA 2 and 3 to that of VLA 3, thus there may be a common rotational motion for all the sources present.

### 4.2.2 3D modelling

There are two best-fitting models (models with the same overall ranking), but the only difference is the outflow reference velocity. Table 8 lists the parameters of the best-fitting models for the CH$_3$CN $J = 12-11, K = 2$ to 5 molecular line emission for a model with a stellar mass of 35 M☉. Unless otherwise stated, we will refer to these models as the best-fitting ones as the stellar mass matches the mass derived from the source luminosity. The best-fitting models have outflow reference velocities $v_0 = 0$ and $1 \text{ km s}^{-1}$. However, this parameter may not be well constrained because the observations do not show any clear evidence of outflow motion in the first moment map, e.g. online Fig. C1. Hence, for simplicity, we plot the results for the model with zero outflow reference velocity. Fig. 6 shows the pv maps of the best-fitting model for each line and cut direction. The model reproduces relatively well the angular extent of the pv maps, but not the spectral extent. The zeroth moment slices in Fig. 14 also show that the angular extent is well fitted.

Fig. 15 shows the spectra at the peak of the zeroth moment map, which is equivalent to the spectra at zero offset in the pv maps, for each K level. The best-fitting models do not match the width and the peak of the observed spectral line, mainly because the model infall and rotation velocities are high enough to produce double-peaked spectra. Fig. 15 shows that models with lower stellar masses are needed to match them, but these models are not physically realistic given the stellar luminosity. In general, the abundance scales the line intensity. The best-fitting models prefer low abundances because larger values start to overestimate the values at the line wings. On the other hand, the observed lines are skewed towards higher velocities whilst the model lines are more symmetric.

The best-fitting models have a turbulent width of $2.0 \text{ km s}^{-1}$. Although this is in the upper end of our model grid, the lines are already wide due to the large rotation/infall velocities. The third model in the tally has a turbulent width of $1.5 \text{ km s}^{-1}$.

# 5 DISCUSSION

## 5.1 70-μm morphology

Our detailed radiative transfer modelling has shown that the bipolar outflow cavities play an important role in reproducing the 70-μm observations. The modelling has allowed the fit of the extended emission at both this wavelength and the SED. The observed and model 70-μm images are elongated along the cavity axis. The results in Figs 10 and 12 imply that the warm/hot dust in the envelope cavity walls bring an important contribution to the emission at 70 μm. Furthermore, the temperature distribution in Fig. 16 shows that cavity walls are being heated by the stellar radiation escaping through the outflow cavity more than other regions of the envelope. The figure also shows that the dust in the cavity is hotter than the walls; however, there is less dust in the cavity. Hence the contribution of the envelope cavity wall is more important at 70 μm, which is also shown in the SED.

In the models of Zhang, Tan & McKee (2013) the 70-μm emission is extended along the outflow direction. The origin of the extended emission in their models is due to dust in the outflow at distances larger than the core radius, otherwise, the emission is dominated by the cavity walls of the envelope. They used a wide cavity angle (~100°) thus at an inclination of 25° the observer is looking directly to the source, hence the emission of the cavity walls in their models is relatively symmetrical. In comparison, in our best-fitting model, the line of sight passes through a section of the envelope.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Results for each parameter from the 1D LTE modelling of the NOEMA compact configuration CH<sub>3</sub>CN *J* = 12 – 11 transition observations. The crosses mark the position of the continuum sources in Table 3. The beam is shown in the lower left-hand corner (black ellipse).</figcaption></figure>

**5.2 Rotation of the envelope**

Our analysis of methyl cyanide interferometric observations shows that they are tracing rotating gas at intermediate scales (~5000 au) in the envelope of AFGL 2591. The velocities in the line of sight and the

radiative transfer modelling imply that rotation dominates rather than infall or outflow motions. The direction of the rotation is consistent with the observations of smaller scales made by Wang et al. (2012). They also found that the linear velocity gradient in the HDO line pv
<figure><img src="image_16.png" /></figure>

**Figure 14.** Observed (black line) and models zeroth moment map slices through the peak along the vertical (upper panel) and horizontal (lower panel) axes for *K* = 3. The dashed green line corresponds to the best-fitting model with zero outflow reference velocity. The dotted red line corresponds to models with the same parameters as the best-fitting one but with the temperature step at *T*₀ = 200 K. The dash-dotted blue line shows a model with the same parameters as the best-fitting model but with a constant abundance.

map is 16.2 km s⁻¹ arcsec⁻¹ at a PA = 0°. Similarly, we calculated a linear velocity gradient from the pv maps at PA = 6° for each *K* line and obtained an average value of 13.3 ± 0.5 km s⁻¹ arcsec⁻¹, which is consistent with slower rotation of the envelope in comparison to the disc. Thus the small-scale rotating motion observed in other molecules can be connected to the large-scale motion traced by methyl cyanide.

We estimated a stellar mass of 35 M☉ from the source luminosity. This mass is high enough for the models to produce double-peaked lines which do not fit the observations. The width of the line in these models is dominated by the large velocity gradient produced by a larger stellar mass rather than the non-thermal velocity width. We explored models with lower masses, and determined that for the Ulrich envelope with Keplerian disc model to fit the lines, stellar masses below 15 M☉ are needed (see Fig. 15) which is unrealistic given the luminosity of the source (under this simple scenario where a single protostar would form). Hence the rotation and/or infall motions of the inner envelope are slower than the ones predicted from the Ulrich envelope with a Keplerian disc. A binary or multiple system would provide less luminosity per total stellar mass, which would not fit SED. The slower rotation may be indicative of material being funnelled through the envelope to the disc and/or the result of magnetic braking. Models at lower inclination angles (face-on), i.e. dominated by infalling motions rather than rotation, are discarded by the dust continuum models.

<figure><img src="image_17.png" /></figure>

**Figure 15.** Observed (black line) and best-fitting model peak spectra (dashed green line) for (a) *K* = 2, (b) *K* = 3, (c) *K* = 4, and (d) *K* = 5. The dash-dotted red lines correspond to the best-fitting model with *M*<sub>⋆</sub> = 25 M☉. The dotted blue lines correspond to the best-fitting model with *M*<sub>⋆</sub> = 15 M☉.

<figure><img src="image_18.png" /></figure>

**Figure 16.** Temperature distribution in the mid-plane of the inner 10 000 au region as obtained for the best-fitting model.
Benz et al. (2007) noted that their observations of HCN molecular line emission at the continuum peak were skewed towards the red as in the observations presented here for CH3CN, whilst other lines (e.g. SO) were skewed towards blue. They interpreted the blue-skewed lines as produced by material in the inner region absorbing the red-shifted emission. Their observations, however, trace a smaller region close to the star, thus this can be explained by an optically thick disc. On the other hand, the CH3CN radiative transfer modelling lines do not show any asymmetry close to the line centre. Although our radiative transfer modelling seems to favour models without outflows, this is not well constrained since the observations do not show a clear evidence of tracing outflow material. The width of the lines seems to be different in the outflow cavity than in the envelope. The wider lines are located in the outflow cavity and can be explained as turbulent motions of gas in the outflow. The emission from the region where the red-shifted outflow is located does not show the same trend. It can be argued that the emission is optically thick towards the red-shifted outflow, which is consistent with the inclination of the outflow axis with respect to the line of sight being closer to face-on.

## 5.3 Circumstellar matter distribution

The physical properties derived from our modelling are consistent with the current picture of AFGL 2591, but its density distribution is better constrained. Wang et al. (2012) estimated a disc dust mass between 1 and 3 M<sub>☉</sub> from the total flux at 1.4 mm from PdBI observations, and a disc radius of ~1300 au (at 3.3 kpc) from a Gaussian fit to the 1.4-mm visibilities. We have shown that to fit the 1.3-mm extended configuration visibilities (Fig. 2b) a more massive and bigger disc is necessary. However, such a disc does not match the spatial extent of the 1.3-mm emission well (Fig. 3b). This can be explained by a closer to edge-on disc, which would look more compact in the outflow direction where the extent of the emission differ the most. Such a closer to edge-on disc can be the result of precession changing the PA and inclination angle at the smaller scales, but as argued in the previous section, the inner envelope favours an inclination closer to face-on. Johnston et al. (2013) obtained disc radii, which were also fixed to the centrifugal radius, of 3.5 × 10<sup>4</sup> and 2.7 × 10<sup>3</sup> au for their envelope with disc and control models, respectively. They also could not constrain the disc mass and radius as they did not fit the mm interferometric data.

The inclination with respect to the line of sight is consistent with the values in the literature (30°–45°, e.g. Hasegawa & Mitchell 1995; van der Tak et al. 2006), which are generally found based on geometrical considerations, and it is constrained between 25°–35°. Larger values like those found by Johnston et al. (2013) are ruled out by the modelling since at 40° the red-shifted outflow cavity starts to emit in the near-IR. This discrepancy may be explained by the data and methods used. Although the red-shifted cavity is observed in the models at higher inclinations, the intensity is much lower than the peak of the blue-shifted cavity. Hence the red-shifted cavity emission will be negligible due to the lower surface brightness sensitivity of the 2MASS observations. Johnston et al. (2013) fitted the SED and 2MASS image profiles, thus their opening angle is not well constrained (cf. their fig. 4). Another reason the models rule out higher inclinations is that they require a larger envelope mass to compensate for the lower flux resulting from the higher inclination, and higher inclinations produce a deeper silicate feature which does not fit the mid-IR spectrum.

The inner radius is ~three times smaller than the one found by Preibisch et al. (2003) and Johnston et al. (2013). Although

the sublimation radius was allowed to change in order to obtain a maximum dust temperature of 1600 K in the inner radius, the sublimation radii of the models did not increase noticeably. Smaller inner radii, near the sublimation radius, are also not strongly ruled out by Johnston et al. (2013). In order to analyse the effect of a larger inner radius on the near-IR observations, a modified version of the best-fitting model was calculated with an inner radius twice as large. This model did not improve the fit to the speckle visibilities (χ<sup>2</sup> = 49). In addition, the fit to the JHK-bands observations is worse than the best-fitting model (reduced χ<sup>2</sup> differences of about 46, 53 and 99, respectively).

Fig. 16 shows that the peak temperature at 1000 au scales reaches values between 200 and 300 K inside the outflow cavities, whilst in the denser cavity walls the peak temperature decreases to 200–250 K. Palau et al. (2014) obtained a temperature of 250 K at 1000 au scales from their simultaneous modelling of the SED for λ > 70 µm and SCUBA 450-µm and 850-µm radial profiles. Gieser et al. (2019) obtained a temperature distribution from 1D LTE fitting of CH3CN and H2CO transitions of the combined NOEMA/CORE observations. From their equations, we obtain a value of ~220 K at 1000 au from the 1D LTE modelling. All these values are consistent with the temperature distribution from our dust emission modelling.

The parabolic shape of the cavity can help to explain the wider cavity opening angle close to the star as observed in maser emission (Sanna et al. 2012). The density in the cavity is also consistent with the one expected for an outflow (e.g. Matzner & McKee 1999) and the opening angle is consistent with previous results (e.g. CS line modelling by van der Tak et al. 1999, maser emission from Sanna et al. 2012). Due to the lack of data, the previous modelling of MYSOs used a constant small density (or empty) in the cavity (e.g. de Wit et al. 2010). This is not preferred in the modelling, as a constant density tends to produce more extended emission in the near-IR than observed. For a fixed cavity density distribution exponent, the density is determined mainly by the near-IR observations and the near-/mid-IR regime of the SED. If the gas in the cavity is expanding isotropically, the outflow rate Ṁ<sub>out</sub> can be described as

$$ \dot{M}_{\text{out}} = 4\pi r^2 v_{\text{out}} \rho_{\text{cav}}, \quad (4) $$

with v<sub>out</sub> is the outflow velocity and ρ<sub>cav</sub> is the cavity mass density at radius r. The outflow rate is roughly 10 per cent the accretion rate Ṁ<sub>acc</sub> (e.g. Richer et al. 2000; Kölligan & Kuiper 2018). If the accretion rate is Ṁ<sub>acc</sub> = Ṁ<sub>env</sub>, an outflow velocity v<sub>out</sub> = 260 km s<sup>-1</sup> is obtained by solving equation (4) with the values in Table 8. This velocity is within the range assumed by Preibisch et al. (2003) for the expansion of the loops (100–500 km s<sup>-1</sup>) based on the <sup>12</sup>CO line wings in the observations of van der Tak et al. (1999). For the radiative transfer model cavity density, the outflow velocity should be in the range 200–320 km s<sup>-1</sup> for infall rates within 20 per cent of the best-fitting model.

In general, the 450- and 850-µm images from the modelling are not elongated along the cavity axis. For the best-fitting model, the major-to-minor axial ratio of 2D Gaussians fitted to the model images is ~1 at both wavelengths, instead of ~1.3 as observed (cf. Section 2.3.2). Elongation along the cavity axis was predicted by de Wit et al. (2010) for 350-µm observations of the MYSO W33A, but this is not observed in the modelling presented here at pc scales. However, their source has a larger inclination angle (60° for their best-fitting model) and their cavity has a constant density distribution of ρ ≈ 2 × 10<sup>-20</sup> g cm<sup>-3</sup>. Our results indicate that the elongation may be the result from earlier stages in the core formation and/or later evolution rather than produced by the cavities. To a lesser extent and as stated
in Section 2.3.2, changes in the PSF may also help to explain why the models are not elongated.

Similarly to other molecules, the methyl cyanide abundance needed to be defined as a piecewise function in order to fit the extension of the emission and the ratios between different *K* lines. Fig. 14 shows that a constant abundance model produces much more extended emission. Moving the temperature step from $T_0 = 100$ to 200 K underestimates the peak zeroth moment, thus the *T* = 100–200 K range is the most relevant one in order to fit the peak of the zeroth moment emission. The inclusion of a step at 230–300 K was not needed during the line radiative transfer modelling. This jump should increase the abundance of *N*-bearing molecules in the inner hotter regions (e.g. Rodgers & Charnley 2001). Such temperatures are reached in regions much closer to the star, which are not resolved by the methyl cyanide observations presented here. The gas temperature from the chemical modelling of the high-resolution NOEMA/CORE A+B+D data by Gieser et al. (2019) also does not go above ~220 K at the peak of the emission. Further refining of the radiative transfer grid is also needed in order to further explore the changes in abundance.

## 5.4 Multiwavelength modelling

Here, we discuss several aspects of the modelling procedure. First, we discuss the sensitivity of the model parameters to our observations and the selection of the best-fitting model given the large amount of data available for this work. Then we suggest improvements to the modelling with a focus in an observational aspect to constrain the less sensitive parameters and a theoretical aspect to improve the physical model, and thus obtain a better interpretation of the data for this and other MYSOs.

### 5.4.1 Sensitivity of model parameters

In total, there are 15 free parameters defining each density distribution and heating source for model types A–C without including the dust models, the stellar mass and stellar radius. Additionally, there are three free parameters describing the kinematics of the gas. The inner radius depends on other parameters, and a few of them are not expected to be well-constrained by the modelling based on the observations available. These include the parameters defining the disc flaring and density structure (three parameters), and the stellar temperature. The latter converged to different values for the three types of models presented in Johnston et al. (2013), but their inner radius, which was not well constrained, was a multiple of the sublimation radius, which, in turn, depends on the stellar temperature. Therefore, there are ten free parameters, which can be constrained by the continuum observations (see column ‘constrained by’ in Table 8).

We calculated grids of models for parameters that can be well constrained by the observations available: inclination angle, cavity opening angle, envelope infall rate, disc mass and disc/centrifugal radius. These grids allow us to estimate the uncertainties of the varied parameters. Fig. 17 shows $\chi^2$ maps as a function of tentatively correlated parameters for specific observations.

Fig. 17(a) shows the $\chi^2$ values for 1.3-mm extended configuration visibilities in a grid of models with varying disc mass and radius, and the other parameters equal to the best-fitting model ones. From the figure, we estimate an error of half the parameter resolution in this region of the grid, i.e. errors of 200 au in disc radius and 1 M⊙ in disc mass. We explored the correlation between these parameters and other observations, and found that in general IR observations

(including 70 µm) and the SED prefer smaller disc radii. On the other hand, sub(mm) observations prefer models with bigger discs and higher masses. However, none of these observations constrained these parameters as well as the 1.3-mm extended configuration visibilities.

The inclination angle and cavity opening angle plane in Fig. 17(b) shows that the minima of the SED and *H*-band observations coincide, whilst the *J* and *K* bands also coincide but at higher angles. However, the contour levels show that the latter have flatter $\chi^2$ distributions than the former, and Fig. 7 shows that good fits are obtained for the *J*- and *K*-bands at the best-fitting angles. The error in this case is dominated by the SED, because its $\chi^2$ distribution is more steeper than the other observations as pointed by its contour. We estimate an error of 5° for both angles. For models with constant envelope infall rates, the SED is best fitted by models with higher values for both angles as the luminosity increases. A similar behaviour is observed in the best fits of the *JHK*-bands, but due to the Monte Carlo noise the positions of the minima are more erratic. On the other hand, if the luminosity remains fixed the inclination angle of the best-fitting model to the SED decreases as the envelope infall rate increases, and there is not a clear trend in the values of the cavity opening angle. This is consistent with a deeper silicate feature at higher envelope infall rates (higher envelope masses) which is compensated by a lower inclination angle. The same is observed in the minima for the *JHK*-bands observations.

We did not use the 70 µm to constrain the inclination and cavity opening angles. The 70-µm distribution in Fig. 17(c) shows that this observation is fitted by larger values for both angles. This is probably pushed by the higher fluxes in the blue-shifted cavity as observed in the horizontal slice in Fig. 1. The inclination angle affects mainly the extension of the emission along the outflow cavity whilst the change in opening angle affects both directions but the change is more subtle. However, the increase in these angles does not increase the emission closer to the source as in the horizontal slice. In addition, the density distribution cannot account for e.g. inhomogeneities in the cavity, which may improve the fit. Submm observations are best fitted by models with lower opening angles, except for the 1.3-mm extended configuration visibilities which prefer slightly larger values. The latter are best fitted by slightly larger inclination angles (mainly in the 30°–40° range) than the best-fitting model one. Even though we did not calculate a new grid of disc mass and radius with the best-fitting values for the inclination and cavity opening angles, Fig. 2 shows that a relatively good fit is obtained.

In Fig. 17(d), we show a cut in the stellar luminosity and envelope infall rate for the best-fitting inclination and cavity opening angles. In the figure, the position of the best-fitting model does not coincide with any of the individual observations, but it is in the intersection of the contours of the observations. Its position is at the same distance to the SED and the *H*- and *K*-band observations, hence it provides relatively good fits to them. The distribution of $\chi^2$ values for the *J*-band image is flatter than the other observations. From these grids, we estimate an error in the envelope infall rate of $0.2 \times 10^{-3}$ M⊙ yr$^{-1}$, i.e. roughly 10 per cent of the best-fitting model, and an error of $0.2 \times 10^5$ L⊙ for the stellar luminosity. The minima for the best-fitted models to the 70-µm and 1.3-mm compact configuration data are generally located towards the higher end of luminosity and infall rate in the grid, whilst the *K*-band speckle and the SCUBA radial profiles prefer values in the lower end of the grid of these parameters. The minima of the 1.3-mm extended configuration prefer higher luminosities and lower infall rates than the best-fitting model.

Fig. 18 shows the position of the minima for grids with changing cavity opening angles and inclinations. For a fixed cavity opening
<figure><img src="image_19.png" /><figcaption>**Figure 17.** χ² maps in parameter space. (a) χ² of 1.3-mm extended configuration interferometric visibilities in the disc mass/radius plane. The red contour corresponds to a 50σχ² level with σχ² = 10. The dot with red border marks the minimum value. (b) Minimum χ² value (dots) and 5σχ² (see online Table A4) contour level for the SED (black), and *J* (red), *H* (green), and *K* (blue) observations in the cavity half-opening/inclination angle plane. (c) Herschel 70-μm χ² map in the half-cavity opening/inclination angle plane. The dot with red border marks the minimum value and red contour corresponds to the 5σχ² level with σχ² = 1. (d) Same as (b) but for the luminosity/infall rate plane. The black cross marks the values adopted for the best-fitting model.</figcaption></figure>

angle (Fig. 18a), in the best-fitting models of the SED the luminosity increases and the infall rate decreases as the inclination angle increases. This is expected because as the inclination angle increases the silicate feature is increasingly deeper, whilst a larger opening angle overestimates the fluxes between ~2 and 8 µm because this range is the most affected by the cavity emission (see Fig. 10). Hence, the envelope infall rate must decrease and the luminosity increase to compensate for the loss of flux in the envelope emission. On the other hand, in the best-fitting models of the JHK-band observations the luminosity and the infall rate decrease as the inclination increases. For a fixed inclination angle (Fig. 18b), in the best-fitting models of the SED and the JHK-band observations the infall rate increases as the opening angle increases, and in the case of the SED the luminosity slowly increases.

As part of the visual inspection step, we explored models with different cavity properties. Between models with cavity shape exponents $b_{cav} = 2$ and 1.5, the SED and the near-IR observations are best fitted by the former. Models with lower cavity shape exponents the minimum χ² has a larger inclination angle by 5° between the shape exponents explored. For larger cavity shape exponents, the line of sight crosses a wider section of the envelope for a fixed inclination and cavity opening angle, hence the difference in minimum inclination angle. Models with a steeper density distribution tend to fit the SED best, reaching a factor of ~2 difference between the χ² of a constant and a $p_{cav} = 2$ density distribution. Similarly, the *J-* and *H*-band observations also favour steeper density distributions, whilst

the K-band observation prefer intermediate exponents. None the less, the difference in χ² is not as sharp as for the SED. A decrease in the cavity reference density of a factor of ~2 or larger increases the χ² of the SED by a similar factor, but the fit is less sensitive to changes of 0.2–0.4 × 10⁻²¹ g cm⁻³. To fit the near-IR observations including the speckle interferometry, different combinations of the cavity parameters can produce noticeable changes as these determine the amount scattered emission. Some combinations produce local minima in χ² values for the UKIRT observation fit, with comparable values between these minima. Given that the emission is not homogeneous, different values may be trying to fit different sections of the cavity. In general, cavity reference densities of 2 or more orders of magnitude smaller than the best-fitting values for the SED are required to fit well the speckle visibilities, with increasingly larger density values for flatter density distributions. All in all, it cannot be ruled out that a different combination of the five parameters determining the aspect of the cavity may produce better results as these are probably degenerate, and we cannot break these degeneracies because the models cannot account for the inhomogeneities observed.

As expected, models with higher submm dust opacities, such as WD01, require a lower accretion rate, i.e. lower envelope mass, to match the submm points, but still within ~20 per cent of the best-fitting model value. Another parameter that determines the envelope mass is its radius, which is constrained mainly by the 850-µm intensity radial profile. Values between ~1.5–2.5 × 10⁵ au fit this profile at least up to 40 arcsec. The simulated observation did not
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Position of the minimum $\chi^2$ in the luminosity–envelope infall rate plane for different (a) inclination angles and (b) cavity opening angles. The colour of the circles represent the SED (black), and *J* (red), *H* (green), and *K* (blue) observations. The sizes of the circles increase with inclination angles (25$^\circ$, 35$^\circ$, 45$^\circ$) and cavity opening angles (25$^\circ$, 28$^\circ$, 30$^\circ$).</figcaption></figure>

take into consideration chopping during the real observations, which can make the observed radius smaller.

Although models with different flaring, scale heights, and vertical density exponents were calculated (cf. online Appendix A1.3), none of them improved the overall fit. Since most of the observations are not sensitive enough to these parameters, an accurate validity range cannot be given.

Stellar parameters are difficult to constrain since the protostellar emission cannot be observed directly in the near-IR. By taking into account the ratios of dust extinction/emission and scattering from our best-fitting model, we tried to constrain the stellar temperature by using the highest resolution near-IR point source data. However, the continuum is still dominated by photons scattered or emitted by dust, so this did not allow us to constrain the stellar temperature due to the low amount of direct stellar photons predicted by the model.

On the other hand, the line fitting is affected mainly by the abundance and line width for a fixed stellar mass. This was taken into account during the developing of the model grids. We estimate that the data can be sensitive to a finer grid of abundance values, with steps of e.g. 0.1 × 10$^{-8}$. The linewidth can also be better explored by a finer grid, but with lines whose width is not dominated by the gas kinematics.

### 5.4.2 Best-fitting model selection

In our analysis, we have used two methods to select the best-fitting models: one based on the ranking of $\chi^2$ values of each observation to select a seed model from the visual fit to dust continuum observations

(Appendix A1.2, online) and to select the line best-fitting model (Appendix A3, online), and another method where we develop a grid of models to fit specific dust continuum observations and parameters (Appendix A2, online). The former ranking method has been used by a few authors rather than using a e.g. $\chi^2$ goodness of fit statistic. Jørgensen, Schöier & van Dishoeck (2002), who combined SED and submm intensity profiles, argue that the $\chi^2$ values cannot be combined because the observations are not completely independent, thus it is not statistically correct to combine them. They also argue that intensity profiles and SEDs constrain different parameters of the density distribution. Williams, Fuller & Sridharan (2005) add that the $\chi^2$ values are in different numerical scales, thus they should be compared in an ordinal way. It can also be added that e.g. radial profiles can have a number of data points comparable to the SED, thus certain wavelengths will carry more weight in a total $\chi^2$. In order to explore this further, we show in Table 9 the continuum best-fitting models for the visual inspection step selected by using a total reduced $\chi^2$ and an average reduced $\chi^2$.

The best-fitting model as selected by the total reduced $\chi^2$ has a hotter star (36 000 K) with the same luminosity, a larger envelope ($R_{\text{out}} = 270000$ au) and a slight higher accretion rate ($M_{\text{env}} = 1.4 \times 10^{-3} \, \text{M}_{\odot} \, \text{yr}^{-1}$) than the best-fitting model as selected by ranking during the visual inspection step. Their discs are the same except for the flaring and vertical density exponents, and their cavities have the same density distribution. As was shown in Section 4.1 and in Fig. 17 higher inclination and opening angles are preferred by the 70-µm observations but not by the *J*- and *H*-band ones, and the SED. It can also be seen in Table 9 that better $\chi^2$ values than in the best-fitting model by ranking are obtained in those observations with a higher number of points, namely the UKIRT *K*-band and 70-µm observations. Hence, this method for selecting the best-fitting model is dominated by the observations with a larger number of points.

The best-fitting model selected by the average $\chi^2$ is a type C model equal to the best-fitting type C model selected by ranking (see online Table D1), except for the envelope radius (200 000 au). As argued, this type of models do not fit the near-IR observations and the SED. In general, this method gives more weight to observations, which have in general a high $\chi^2$ value, namely the 1.3-mm NOEMA observations. It also reflects the difference in $\chi^2$ scale between the different observations. If all the observations were well fitted ($\chi^2 \sim 1$) then this best-fitting model should be the same as the ranking one.

In general, the method used in the visual inspection fitting is independent of the intrinsic $\chi^2$ scale of each observations and the number of points of each one. Overall the best-fitting model selected by ranking provides a good fit for key observations and its total and average $\chi^2$ values are relatively close to the best-fitting ones selected by these two $\chi^2$. The disadvantage of the ranking is that it reduces the difference between models with a large difference in $\chi^2$ or increases it for models with similar $\chi^2$.

On the other hand, in the grid modelling step of Fig. 11 (see also Appendix A2, online), we determined the best-fitting model from the intersection of confidence regions in parameter space built from $\chi^2$ values for different dust continuum observations. This procedure is commonly used to constrain parameters in cosmological models. The advantage of this method is that it does not depend on the intrinsic scale of the $\chi^2$ values and does not reduce the differences in $\chi^2$ value between models. However, prior knowledge of the effect of each parameter in the fit to each observation is needed, and it requires more computing time to develop a grid. Hence combining this method with a prior filter from our visual inspection fitting improved the $\chi^2$ values of key observations.
<figure><img src="image_21.png" /><figcaption>**Figure 19.** SED of the cm/submm regime as observed at different scales. Blue squares fluxes correspond to single-dish observations and are presented in Table 6. The orange triangle corresponds to the 1.3-mm extended configurations flux in Table 3, and the green triangles come from mm interferometric observations from van der Tak et al. (1999), Johnston et al. (2013), and Wang et al. (2012). Red circles are radio/mm interferometric observations from Campbell (1984), Tofani et al. (1995), Trinidad et al. (2003), van der Tak & Menten (2005), and Johnston et al. (2013). The slopes of the fitted lines are listed in the upper left-hand corner. If the emission is from dust, the dust emissivity slope is $\beta = \alpha - 2$.</figcaption></figure>

### 5.4.3 Dust size distribution and grain growth

In order to explore the change in dust from within the envelope or disc, we derived spectral indices based on dust emission at different scales from single-dish and interferometric observations in the radio/submm regime. The former traces cold dust located in the large-scale envelope, whilst the latter trace scales of a few 1000 au. The dust optical depth is assumed to follow a power law $\tau \propto \nu^{\beta}$ in this regime (optically thin). If it is further assumed that the emission is in the Rayleigh-Jeans regime, the fluxes are then $F_{\nu} \propto \nu^{\beta+2} \propto \nu^{\alpha}$. Fig. 19 shows the observations made at different scales and the indices $\alpha$ of the power law fitted to the single dish, mm, and radio interferometry data individually.

Dust models with shallower emissivity indices (e.g. WD01, see online Table A3) fit the $\lambda > 70$-µm regime of the SED better than dust with steeper spectral indices (e.g. KMH, OHM92), and allow the fit of visibilities at smaller baselines in the 1.3/1.4 mm with smaller discs that better reproduce the 1.4-mm emission. This is supported by the emissivity index derived from single dish observations ($\beta = 1.2$ from Fig. 19). However, shallower emissivity index models, i.e. models with larger grains (cf. online Table A3), are not good at fitting the mid-IR regime and show an excess in emission along the cavity in the near-IR images because larger grains scatter more photons.

The dust emissivity index derived from interferometric observations is $\beta = 0.2$. This shallower emissivity index implies that larger grains are responsible for the mm disc emission, thus it may be

evidence of grain growth in the denser regions, i.e. the disc, of the MYSO. However, contamination from free-free emission can contribute to this index, hence a higher emissivity index cannot be ruled out. Although the best-fitting model from the visual inspection step favours a disc with smaller grains, good fits to the 1.3-mm extended configuration visibilities can be obtained with larger grains (W02-3 grains), albeit reducing the disc mass. On the whole, both types of grains actually do a good job at reproducing the observations, so there is not strong evidence that grain growth is needed to explain the data.

The spectral index of radio/mm observations and 3.6-cm observations (e.g. Trinidad et al. 2003; Johnston et al. 2013) points towards the presence of a jet or wind with $\alpha \sim 0.6$. Hence, the spectral index from mm interferometry may be steeper if the contribution of the jet/wind emission is subtracted, which will depend on the turnover frequency of the jet emission.

The presence of different dust properties to those used here is also supported by modelling of near-IR polarimetric observations of Simpson et al. (2013). They found that elongated dust grains rather than spherical ones can reproduce their observations better. Therefore, the dust properties may be different depending on the physical conditions at different positions, e.g. dust in the cavity walls is being shocked by the outflows or dust in colder regions may develop ice mantles. This may allow dust with $\beta = 1.5$ to be present in colder regions of the envelope, thus improving the fit of the submm observations, and dust with $\beta = 2.0$ to exist in hotter or shocked regions like the cavity walls, which may improve the mid-IR fit.

### 5.4.4 The Ulrich solution limitations

The Ulrich (1976) analytical solution for the accretion problem is able to reproduce some of the observed features. However, we could not find one model that can fit all the observations together. In particular, the results of the continuum and CH$_3$CN line modelling are contradictory. For instance, the stellar mass required by the Ulrich envelope with Keplerian disc velocity model is <15 M$_{\odot}$, which is much smaller than the one inferred from the source luminosity. This problem cannot be solved by changing the inclination angle, as lower inclination angles do not provide good fits to the J- and K-band observations. The slower observed rotation than the one suggested by the theoretical model may require the effects of magnetic braking to slow down the rotation of the disc and inner envelope. These limitations on the physics included in the Ulrich solution do not allow us to obtain better results.

Additionally, the simplicity of the models, where the clumpiness of the cavity or the presence of nearby sources are not taken into account, also limits the quality of the fit. The high-resolution 1.3-mm extended configuration observations in Fig. 3(b) recover more extended emission than the previous 1.4-mm observations from Wang et al. (2012), and show that the emission is more clumpy/less axisymmetric than the previous observations. This has also been observed in other MYSOs, e.g. high-resolution Atacama Large Millimeter/submillimeter Array (ALMA) observations of W33A (Maud et al. 2017). In the case of W33A filaments resembling spiral arms has been interpreted as accretion flows (Izquierdo et al. 2018). The velocity distribution in these accretion flows can be modelled with an Ulrich envelope by selecting some of its streamlines (e.g. Tobin et al. 2012, for low-mass star-forming regions). Hence, a full Ulrich envelope, rather than selecting some streamlines, is less likely.

Simulations with magnetic fields and radiation transfer show promising results. The radiation-magnetohydrodynamics (RMHD)
simulations of early stages in the formation of high-mass stars calculated by Commerçon, Hennebelle & Henning (2011) show that magnetic fields and radiation can reduce the fragmentation during early stages in the formation of massive stars. The collapse of a $10^3 M_\odot$ rotating molecular cloud investigated by Peters et al. (2011) through RMHD simulations also shows that magnetic fields can reduce fragmentation and that they are important in transporting angular momentum at local scales. Furthermore, the RMHD simulations of Myers et al. (2013), which explore a larger time evolution of a collapsing core with an initial power law density distribution index of -1.5, show that magnetic fields can remove angular momentum from the infalling material through magnetic braking whilst still forming a massive star ($20 M_\odot$). Their models can still form an accretion disc but their centrifugal radii are <50 au, thus smaller than the one of AFGL 2591. Larger accretion discs can be obtained when non-ideal MHD effects are taken into account (Kölligan & Kuiper 2018). These RMHD models also show that the matter is accreted into the disc from filaments which have collapsed within the envelope at scales of 5000 au in the mid-plane. However, these models do not include feedback from an outflow.

The non-magnetic simulations by Klassen et al. (2016) produce a Keplerian disc but they focus on 1000 au scales, thus in the transition from the envelope to disc. They also have a higher rotational to kinetic energy ratio than Myers et al. (2013) by a factor of ~5–30, and impose an initial solid body rotation. In the 3D radiation-hydrodynamic simulations of Meyer et al. (2018), Keplerian discs are produced which develop spiral arms and fragment into a multiple stellar system. The 2D simulations of Kuiper et al. (2015, 2016) and Kuiper & Hosokawa (2018) show the effect of radiation on the outflow cavity and its widening from an initially slowly rotating core. The cavity distributions in these models resemble the near-IR observations of AFGL 2591. However, they do not include magnetic fields. Newer simulations are including the effects of an increasing number of physical processes, e.g. in the Kuiper & Hosokawa (2018) simulations the effects of photoionization from a potential HII region are included. These and future simulations have the potential to describe better the physical properties in MYSOs from high-resolution observations.

### 5.4.5 A recipe for future multiwavelength modelling

The diagram in Fig. 20 summarizes the steps necessary for modelling multiwavelength data with 2D radiative transfer models in order to obtain better constraints to the model physical parameters from specific observations. The modelling is performed in three steps with increasing refinement of key parameters at each step and complexity in the modelling, and each step is divided in a series of sub-steps. The sub-steps are in sequential order, although some of them can be swapped or skipped depending on the available data. Here, we discuss the main features and caveats of each step.

The aim of the first step is to obtain initial values and constraints for the 2D model parameters in Table 8. The observations for each parameter in Fig. 20 are in order of preference. Some of these parameters will then be further constrained in the next steps. The initial value for the luminosity can be obtained from model templates (e.g. Mottram et al. 2011a). These model templates can account for rotationally flattened envelopes and cavities (e.g. Robitaille 2017), hence obtaining an initial range of inclination angles and cavity properties can help to constrain the templates. We have listed the disc shape (disc flaring exponent, scale height) in two sub-steps. The first one is for cases of discs closer to edge-on using imaging whilst

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 20.** Steps for modelling and bringing up constraints from MYSOs observations. Parameters under 'cavity shape' are opening angle and cavity shape exponent. Parameters for 'disc shape' are scale height and flaring exponent. The item 'outflows' refers to molecular line observations of outflows. Simple modelling includes fitting with templates or with models such as 2D Gaussians.</figcaption></figure>

the second relies on the modelling of mid-IR visibilities (e.g. di Folco et al. 2009). Although measuring these disc properties, and to a lesser extent the disc radius, from mid-IR observations is desirable as they are tracing optically thick emission from the surface of the disc, some MYSOs may still only be partially resolved in the mid-IR in high-resolution observations (e.g. de Wit et al. 2009). Envelope and disc masses measured from dust emission are obtained by estimating
a dust temperature from the SED. These are a good approximation in terms of order of magnitude for the next step. Some deeply embedded MYSOs may be weak in the near-IR, hence to derive the properties of the outflow cavity shape (inclination angle, opening angle, shape exponent) molecular line emission tracing the outflow (e.g. CO) can be used instead. However, precession of the outflow axis can affect the determination of the cavity geometry, PA, and inclination angle. In such cases, high-resolution mm interferometric observations tracing the disc can provide further constraints. In general, most geometric properties are degenerate due to projection effects, except at extreme inclination angles.

In Step 2, we have integrated the main relations found in our modelling and added further steps to constrain other parameters. All the sub-steps in this step constrain the physical parameters through 2D axisymmetric modelling and are relatively independent from each other. Some sub-steps iterate between parameter-observations relations in order to optimize the modelling, however, *n*-dimensional optimization, with *n* being the number of parameters in the sub-step, can also be performed. Additional parameters constraining the velocity distribution are introduced in this step after modelling the continuum data. Initial values for these parameters can be derived from the luminosity (stellar mass) or analysis of the spectra (outflows).

The objective of the final step is to model 2D axisymmetric features of specific observations in order to shrink the area in $\chi^2$ space in the relations of the previous step. Hence, the sub-steps in the previous step must be repeated with the new model and an additional iteration must be performed to constrain the parameters of the newly introduced features (see example in Fig. 20). Since the sub-steps are relatively independent, the introduction of other 2D structures may not affect all the relations. For instance, introducing entrainment features to better model near-IR observations will likely not affect the envelope radius and disc properties. On the other hand, introducing a torus to explain mid-IR features (e.g. de Wit et al. 2011) can affect the disc properties as well as the envelope mass. Determining the orientation of the 2D features should be straight forward because the orientation angles (PA and inclination) were constrained in the previous steps. This step requires that the comparison between model and observations is performed on a pixel-by-pixel basis for images in order to increase the number of data points needed to fit the increasing number of parameters in the model. Some of the additional model features may be geometrical rather than physical, hence this step may only constrain the physical parameters in the previous one. Non-axisymmetric features can then be introduced to the models to optimize the relationships in Step 2 (see Section 5.4.4 for examples).

# 6 CONCLUSIONS

We used spatially resolved dust continuum images between 1 and 1400 µm and the SED to study the dust density and temperature distributions, and methyl cyanide interferometric observations to study the inner envelope velocity distribution of the prototypical MYSO AFGL 2591. We performed a multistep radiative transfer modelling to constrain the physical model parameters and determined correlations between specific observations and model parameters. We show for the first time that the Herschel 70-µm emission is extended in the outflow direction and explain its extended nature as the result of dust emission from the outflow cavity, whose walls are heated by the radiation escaping through the less dense cavity medium. We found that the inner envelope is rotating and this rotation is slower than the disc.

We have attempted to fit a large number of observations with a unified physical model. The best-fitting model consists of an Ulrich envelope, a flared disc, and bipolar outflow cavities with properties consistent with previous results. The presence of a disc is needed in order to fit NOEMA 1.3-mm extended configuration visibilities. We derive a disc mass of 6 M$_{\odot}$ and 2200 au radius. These values are relatively independent of the parameters determining the envelope/outflow cavities density distribution. The 450-µm and K-band speckle interferometry observations were not fitted well by the models.

This model cannot explain the observed CH$_3$CN pv maps, because the model lines are wider due to the high-velocity gradients of the Ulrich envelope with Keplerian disc model. In addition, the model lines are not skewed like the observed ones. Accretion flows and magnetic braking may be able to explain why the inner envelope rotation is slower than predicted by the Ulrich envelope.

The physical limitations of the Ulrich solution, together with the simplicity of the models, do not allow an improvement over what we have done in this paper. 3D self-consistent theoretical models should be tested in order to improve the fit to the observations. Additional improvement to the modelling and high-resolution multiwavelength observations may also be needed to constrain these parameters. Future papers including higher resolution data from CORE and other surveys will allow a more complete study of MYSOs to test these models.

# ACKNOWLEDGEMENTS

The authors would like to thank T. Preibisch, J. P. Simpson and K.-S. Wang, who kindly provided their data. FO acknowledges the support of CONICYT Becas-Chile 72130407. FO and HRVC acknowledge the support of the Ministry of Science and Technology of Taiwan, project no. 109-2112-M-007-008-. AA, HB, JCM, and CG and SS acknowledge support from the European Research Council under the European Community's Horizon 2020 framework program (2014-2020) via the ERC Con-solidator Grant 'From Cloud to Star Formation (CSF)' (project number 648505). HB also acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) via the Collaborative Research Center (SFB 881) 'The Milky Way System' (subproject B1). AP acknowledges financial support from CONACyT and UNAM-PAPIIT IN113119 grant, México. DS acknowledges support by the Deutsche Forschungsgemeinschaft through SPP 1833: 'Building a Habitable Earth' (SE 1962/6-1). RK 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 Foundation (DFG) under grant no. KU 2849/3-1 and KU 2849/3-2. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration, 2013). Part of this work is based on archival data, software or online services provided by the ASI Science Data Center (ASDC).

# DATA AVAILABILITY

The data underlying this paper will be shared on a reasonable request to the corresponding author. Parts of these data were derived from: WFCAM Science Archive (http://wsa.roe.ac.uk/), Herschel Science Archive (http://archives.esac.esa.int/hsa/whsa/), SCUBA Legacy Catalogues (http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/community/scubalegacy/). The K-band speckle interferometric data were provided by T. Preibisch by permission.
REFERENCES

Ahmadi A. et al., 2018, A&A, 618, A46

Benz A. O., Stäuber P., Bourke T. L., van der Tak F. F. S., van Dishoeck E. F., Jørgensen J. K., 2007, A&A, 475, 549

Beuther H. et al., 2018, A&A, 617, A100

Boonman A. M. S., Stark R., van der Tak F. F. S., van Dishoeck E. F., van der Wal P. B., Schäfer F., de Lange G., Laauwen W. M., 2001, ApJ, 553, L63

Bruderer S., Benz A. O., Doty S. D., van Dishoeck E. F., Bourke T. L., 2009, ApJ, 700, 872

Bruderer S., Benz A. O., Stäuber P., Doty S. D., 2010, ApJ, 720, 1432

Campbell B., 1984, ApJ, 287, 334

Casali M. et al., 2007, A&A, 467, 777

Commerçon B., Hennebelle P., Henning T., 2011, ApJ, 742, L9

Cunningham A. J., Klein R. I., Krumholz M. R., McKee C. F., 2011, ApJ, 740, 107

Davies B., Hoare M. G., Lumsden S. L., Hosokawa T., Oudmaijer R. D., Urquhart J. S., Mottram J. C., Stead J., 2011, MNRAS, 416, 972

de Wit W. J. et al., 2009, A&A, 494, 157

de Wit W. J., Hoare M. G., Oudmaijer R. D., Lumsden S. L., 2010, A&A, 515, A45

de Wit W. J., Hoare M. G., Oudmaijer R. D., Nürnberger D. E. A., Wheelwright H. E., Lumsden S. L., 2011, A&A, 526, L5

di Folco E., Dutrey A., Chesneau O., Wolf S., Schegerer A., Leinert C., Lopez B., 2009, A&A, 500, 1065

Di Francesco J., Johnstone D., Kirk H., MacKenzie T., Ledwosinska E., 2008, ApJS, 175, 277

Duarte-Cabral A., Bontemps S., Motte F., Hennemann M., Schneider N., André P., 2013, A&A, 558, A125

Gieser C. et al., 2019, A&A, 631, A142

Griffin M. J. et al., 2010, A&A, 518, L3

Haemmerlé L., Peters T., 2016, MNRAS, 458, 3299

Hasegawa T. I., Mitchell G. F., 1995, ApJ, 451, 225

Hatchell J., Fuller G. A., Millar T. J., Thompson M. A., Macdonald G. H., 2000, A&A, 357, 637

Hoare M. G., 2002, in Crowther P., ed., ASP Conf. Ser. Vol. 267, Hot Star Workshop III: The Earliest Phases of Massive Star Birth. Astron. Soc. Pac., San Francisco, p. 137

Hoare M. G., Franco J., 2007, Ap&SS, 1, 61

Hosokawa T., Yorke H. W., Omukai K., 2010, ApJ, 721, 478

Izquierdo A. F., Galván-Madrid R., Maud L. T., Hoare M. G., Johnston K. G., Keto E. R., Zhang Q., de Wit W.-J., 2018, MNRAS, 478, 2505

Johnston K. G., Shepherd D. S., Robitaille T. P., Wood K., 2013, A&A, 551, A43

Jørgensen J. K., Schöier F. L., van Dishoeck E. F., 2002, A&A, 389, 908

Kaźmierczak-Barthel M., van der Tak F. F. S., Helmich F. P., Chavarría L., Wang K.-S., Ceccarelli C., 2014, A&A, 567, A53

Kaźmierczak-Barthel M., Semenov D. A., van der Tak F. F. S., Chavarría L., van der Wiel M. H. D., 2015, A&A, 574, A71

Keto E., Rybicki G., 2010, ApJ, 716, 1315

Klassen M., Pudritz R. E., Kuiper R., Peters T., Banerjee R., 2016, ApJ, 823, 28

Kölligan A., Kuiper R., 2018, A&A, 620, A182

Krumholz M. R., McKee C. F., Klein R. I., 2005, ApJ, 618, L33

Krumholz M. R., Klein R. I., McKee C. F., 2012, ApJ, 754, 71

Kuiper R., Hosokawa T., 2018, A&A, 616, A101

Kuiper R., Yorke H. W., 2013, ApJ, 772, 61

Kuiper R., Klahr H., Beuther H., Henning T., 2010, ApJ, 722, 1556

Kuiper R., Yorke H. W., Turner N. J., 2015, ApJ, 800, 86

Kuiper R., Turner N. J., Yorke H. W., 2016, ApJ, 832, 40

Lada C. J., Thronson H. A., Jr, Smith H. A., Schwartz P. R., Glaccum W., 1984, ApJ, 286, 302

Lucas P. W. et al., 2008, MNRAS, 391, 136

Lumsden S. L., Hoare M. G., Urquhart J. S., Oudmaijer R. D., Davies B., Mottram J. C., Cooper H. D. B., Moore T. J. T., 2013, ApJS, 208, 11

Lutz D., 2012, Herschel Internal Document PICC-ME-TN-033, PACS Photometer Point Spread Function. available at https://www.cosmos.esa.int/documents/12133/996891/PACS+photometer+point+spread+function

Matzner C. D., McKee C. F., 1999, ApJ, 526, L109

Maud L. T., Hoare M. G., Gibb A. G., Shepherd D., Indebetouw R., 2013, MNRAS, 428, 609

Maud L. T., Hoare M. G., Galván-Madrid R., Zhang Q., de Wit W. J., Keto E., Johnston K. G., Pineda J. E., 2017, MNRAS, 467, L120

McDonald M., Courteau S., Tully R. B., Roediger J., 2011, MNRAS, 414, 2055

Meyer D. M. A., Kuiper R., Kley W., Johnston K. G., Vorobyov E., 2018, MNRAS, 473, 3615

Minchin N. R., Hough J. H., McCall A., Aspin C., Hayashi S. S., Yamashita T., Burton M. G., 1991, MNRAS, 251, 508

Molinari S. et al., 2010, PASP, 122, 314

Morales E. F. E., Mardones D., Garay G., Brooks K. J., Pineda J. E., 2009, ApJ, 698, 488

Motte F., Bontemps S., Schilke P., Schneider N., Menten K. M., Broguière D., 2007, A&A, 476, 1243

Motte F. et al., 2010, A&A, 518, L77

Mottram J. C. et al., 2011a, A&A, 525, A149

Mottram J. C. et al., 2011b, ApJ, 730, L33

Myers A. T., McKee C. F., Cunningham A. J., Klein R. I., Krumholz M. R., 2013, ApJ, 766, 97

Olguin F. A., Hoare M. G., Wheelwright H. E., Clay S. J., de Wit W.-J., Rafiq I., Pezzuto S., Molinari S., 2015, MNRAS, 449, 2784

Palau A. et al., 2014, ApJ, 785, 42

Palau A., Sánchez Contreras C., Sahai R., Sánchez-Monge Á., Rizzo J. R., 2013, MNRAS, 428, 1537

Parkin E. R., Pittard J. M., Hoare M. G., Wright N. J., Drake J. J., 2009, MNRAS, 400, 629

Peters T., Klessen R. S., Mac Low M.-M., Banerjee R., 2010, ApJ, 725, 134

Peters T., Banerjee R., Klessen R. S., Mac Low M.-M., 2011, ApJ, 729, 72

Pilbratt G. L. et al., 2010, A&A, 518, L1

Poetzel R., Mundt R., Ray T. P., 1992, A&A, 262, 229

Poglitsch A. et al., 2010, A&A, 518, L2

Preibisch T., Balega Y. Y., Schertl D., Weigelt G., 2003, A&A, 412, 735

Richer J. S., Shepherd D. S., Cabrit S., Bachiller R., Churchwell E., 2000, in Mannings V., Boss A. P., Russell S. S., eds, Protostars and Planets IV. Univ. of Arizona Press, Tucson, AZ, p. 867

Robitaille T. P., 2011, A&A, 536, A79

Robitaille T. P., 2017, A&A, 600, A11

Rodgers S. D., Charnley S. B., 2001, ApJ, 546, 324

Roy A. et al., 2011, ApJ, 727, 114

Roy A. et al., 2014, A&A, 562, A138

Russeil D., Zavagno A., Motte F., Schneider N., Bontemps S., Walsh A. J., 2010, A&A, 515, A55

Rygl K. L. J. et al., 2012, A&A, 539, A79

Sadavoy S. I. et al., 2013, ApJ, 767, 126

Sanna A., Reid M. J., Carrasco-González C., Menten K. M., Brunthaler A., Moscadelli L., Rygl K. L. J., 2012, ApJ, 745, 191

Schneider N. et al., 2016, A&A, 591, A40

Simpson J. P. et al., 2013, MNRAS, 435, 3419

Smith R. J., Longmore S., Bonnell I., 2009, MNRAS, 400, 1775

Stäuber P., Doty S. D., van Dishoeck E. F., Benz A. O., 2005, A&A, 440, 949

Tamura M., Yamashita T., 1992, ApJ, 391, 710

Tobin J. J., Hartmann L., Bergin E., Chiang H.-F., Looney L. W., Chandler C. J., Maret S., Heitsch F., 2012, ApJ, 748, 16

Tofani G., Felli M., Taylor G. B., Hunter T. R., 1995, A&AS, 112, 299

Trinidad M. A. et al., 2003, ApJ, 589, 386

Trinidad M. A. et al., 2013, MNRAS, 430, 1309

Ulrich R. K., 1976, ApJ, 210, 377

van der Tak F. F. S., Menten K. M., 2005, A&A, 437, 947

van der Tak F. F. S., van Dishoeck E. F., Evans II N. J., Bakker E. J., Blake G. A., 1999, ApJ, 522, 991

van der Tak F. F. S., Walmsley C. M., Herpin F., Ceccarelli C., 2006, A&A, 447, 1011

van der Tak F. F. S., Black J. H., Schöier F. L., Jansen D. J., van Dishoeck E. F., 2007, A&A, 468, 627
van Dishoeck E. F., Blake G. A., 1998, ARA&A, 36, 317

Vázquez-Semadeni E., Palau A., Ballesteros-Paredes J., Gómez G. C., Zamora-Avilés M., 2019, MNRAS, 490, 3601

Wang K.-S., van der Tak F. F. S., Hogerheijde M. R., 2012, A&A, 543, A22

Williams S. J., Fuller G. A., Sridharan T. K., 2005, A&A, 434, 257

Zhang Y., Tan J. C., McKee C. F., 2013, ApJ, 766, 86

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

**Figure A1.** Peak spectrum of the CH<sub>3</sub>CN *J* = 12 – 11 observations showing the *K* transitions.

**Figure C1.** NOEMA compact configuration first moment map for the observed CH<sub>3</sub>CN *J* = 12 – 11 *K* = 3 line (a) and from the overall best-fitting model (b).

**Table A1.** Density distributions.

**Table A2.** Parameter ranges used in the visual inspection and grid fitting steps.

**Table A3.** Dust libraries.

**Table A4.** Standard deviation of χ² values from a random model run five times.

**Table A5.** Values of the line modelling parameter grid.

**Table D1.** Parameters of the continuum best-fitting models by density distribution type.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

This paper has been typeset from a TeX/Latex file prepared by the author.