© A. G. Barr et al. 2022

1 Introduction

Circumstellar disks play a key role in the star formation process (Nakano et al. 1995; Jijina & Adams 1996; Bonnell & Bates 2006) as angular momentum conservation of a collapsing cloud leads in a natural way to a disk. Subsequent viscous dissipation and a disk wind, slowly redistribute angular momentum, resulting in accretion of disk material onto the protostar. Massive protostars evolve through several stages - IR dark clouds, Hot Cores, hyper compact HII regions (HCHII), and ultracompact H II regions (UCHII; Zinnecker & Yorke 2007; Beuther et al. 2007; Tan 2017). The Hot Core phase has been associated with the disk accretion phase (Johnston et al. 2015; Maud et al. 2019; Barr et al. 2020), although to what extent this is the case remains to be established.

Ro-vibrational absorption lines at infrared wavelengths trace a range of different energies and are particularly sensitive to the warm material close to protostars. The study of these transitions are therefore a very useful tool for probing the physical and chemical conditions of massive star formation and tracing massive young stellar object (MYSO) structure on small scales. Such a tool can help to understand hot cores and the presence of potential disks, planets, outflow or winds, as well as the accretion processes of massive protostars (Barr et al. 2020).

The disk structure and evolution in low and intermediate mass stars is well understood compared to high mass stars, and significant progress continues to be made (Williams & Cieza 2011; Andrews 2020). Additional processes, such as photoevaporation from high radiative forces, make things more complicated in the study of disks around high mass stars (Zinnecker & Yorke 2007). It is expected that in these objects, disks should also play a fundamental role in the build-up of the stellar mass, as in the case for low mass stars, because otherwise the high radiation pressure would halt accretion before a massive star could be formed (McKee & Ostriker 2007). Furthermore, strong evidence for the presence of disks around massive stars lies in the observation of large bipolar outflows (Shepherd & Churchwell 1996a,b; Wu et al. 2005).

Cirumstellar disks are an important step in the progression from star formation to planet formation, therefore understanding them is vital to gain insight into how planets form, and of what material they are made of. Near-IR observations of low mass stars have been shown to be sensitive to the inner disk on scales of a few AU (Najita et al. 2003; Blake & Boogert 2004; Carr & Najita 2011; Mandell et al. 2012; Gibb & Horne 2013), providing insight to the physical and chemical conditions during planet formation. Planets have yet to be detected around O-type stars with the most massive host stars observed with a planet being 3–5 M_⊙, and less discoveries as the stellar mass increases (Reffert et al. 2015). This likely reflects the short Kelvin-Helmholtz timescales, rapid disk dispersal and high radiation fields of massive stars compared to low mass stars, thus planets have much less time to form (Veras et al. 2020).

At infrared wavelengths disks of massive stars have mainly been studied through atomic line and CO overtone bandhead emission which trace the inner disk (Bik et al. 2006; Davies et al. 2010; Cooper et al. 2013; Fedriani et al. 2020). In a survey of massive protostars, Ilee et al. (2013) modelled CO bandhead emission from 20 massive young stellar objects (YSOs) and find temperatures and densities consistent with emission from an accretion disk, close to the dust sublimation radius. The infrared spectrum of disks around TTauri and Herbig stars is well studied compared to massive YSOs. CO has been extensively studied for both TTauri and Herbig disks (Najita et al. 2003; Blake & Boogert 2004; Salyk et al. 2009; Bast et al. 2011), with H2O and OH also detected in TTauri disks (Pontoppidan et al. 2010; Carr & Najita 2011; Salyk et al. 2011; Banzatti et al. 2017). Herbig disks have been found to consistently exhibit a very low detection rate of H2O compared to TTauri disks, although OH is more readily detected (Pontoppidan et al. 2010; Salyk et al. 2011; Fedele et al. 2012, 2013; Banzatti et al. 2017). Mandell et al. (2012) detected several simple species in emission towards three TTauri disks in the L band with CRIRES and NIRSPEC such as OH, HCN, H2O and C₂H2, also giving upper limits on CH₄ and NH₃. Simple organics have also been detected in emission in the mid-infrared (mid-IR) spectrum of TTauri and Herbig disks with Spitzer (Carr & Najita 2008, 2011; Pascucci et al. 2009; Pontoppidan et al. 2010; Salyk et al. 2011).

Observations of massive stars imply that the mid-IR continuum originates from the dust photosphere of the disk at distances from the star between 50 and 100 AU (Preibisch et al. 2003; Monnier et al. 2009; de Wit et al. 2011; Boley et al. 2013; Frost et al. 2021). However, at mid-IR wavelengths, water and simple organics are seen in absorption towards massive protostars (Evans et al. 1991; Knez et al. 2009; Barentine & Lacy 2012; Rangwala et al. 2018; Dungee et al. 2018; Barr et al. 2018, 2020). Barr et al. (2020) developed a model for AFGL 2591 and AFGL 2136 in which the absorption lines trace gas in an internally heated circumstellar disk at radii of 100-200 AU, with a temperature structure that decreases with scale height, resulting in the presence of absorption lines. We detected HCN and C₂H2 in emission in the L band spectrum of the massive protostar AFGL 2591 and proposed it was due to scattering of continuum photons in the upper layers of the disk photosphere.

In this article we present a follow-up survey of the L band of massive protostars to investigate how widespread emission lines of organics are at these wavelengths, and to understand how they relate to the disk as probed by 4–13 µm absorption lines (Barr et al. 2020). We have selected a wide range of objects in different evolutionary stages: hot cores, ultracompact H II regions and protostars, with particular focus on MonR2 IRS 3 and AFGL 2136. MonR2 IRS 3 differs from AFGL 2591 mainly in the fact that it is a binary, and that typical temperatures are relatively cooler, being on the order of a few hundred K in MonR2 IRS 3 compared to 600–700 K in AFGL 2591. Also the luminosity of MonR2 IRS 3 is an order of magnitude lower than that of AFGL 2591, and the mass is around a factor of 4 lower. AFGL 2136 is very similar to AFGL 2591 in mass, luminosity and temperature of detected molecules, although the inclination angle is 40° (Maud et al. 2019) compared to AFGL 2591 which is face on (van der Tak et al. 1999).

2 Observations and data reduction

All objects were observed as part of program 2021A095 with iSHELL (Rayner et al. 2016) at the Infrared Telescope Facility (IRTF) on Mauna kea. Observations were taken in the L2 mode, at a resolving power of 50000 (6 km s−1), using a slit with of 0.75″. Targets were nodded along the slit to allow for subtraction of sky and telescope emission. The extracted A and B spectra were then averaged. The continuum signal-to-noise (S/N) values per resolution element were >300 for all sources except AFGL 2136 for which only a S/N of 100 was achieved. Flat fields for both science target and standard stars were taken using iSHELL’s internal lamp. Details of the observations are summarised in Table 1.

MonR2 IRS 3 was observed on UT 2021 February 14 from 05:50 to 06:08. Seeing was on average 0.7″ and the source was observed at transit with airmass ranging from 1.13 to 1.15. HR1931 was observed as a standard star immediately before MonR2 IRS 3. HR1931, also known as σ-Ori AB is a O9.5V star with V band magnitude 3.8. It is the brightest star in the σ-Orionis cluster, with three companions. The 0.5–2 µm spectra energy distribution (SED) shows a number of spectral features due to H and He seen in absorption (Simón-Díaz et al. 2015).

MonR2 IRS 2 was observed on UT 2021 February 14 and 15 from 06:19 to 06:58 and 05:37 to 06:26, respectively. HR1948 was used as a standard star on the 14th and HR1931 was used on the 15th. Due to the weakness of MonR2 IRS 2, the achieved S/N was lower than MonR2 IRS 3, and a S/N of 300 was only obtained for part of the spectrum, with most of it lower. HR1948, or ζ Ori Aa, is a O9I star with V band magnitude of 1.9. It is a magnetic hot supergiant and shows evidence for emission lines in the 0.4–4.2 µm SED (Touhami et al. 2010).

We observed Orion BN on February 16 at two different slit angles, from 06:51 to 07:03 and 07:13 to 07:28 at 210° and 235° respectively. HR1948 was used as a standard star on the same night.

AFGL 2136 was observed on UT 2021 July 15 from 07:06:09 to 07:46:56 and 08:15:07 to 08:31:27, with airmass ranging from 1.4 to 1.2. HR7236 was observed in between these two acquisitions as a standard star. HR7236 (λ Aquilae) is a B9V star with a V band magnitude of 3.43. This star has a strong hydrogen Pfund β(7–5) absorption line at 4.654 µm (Bast et al. 2011); we did not detect any evidence for stellar absorption or emission lines for HR7236 in our spectra. Our targeted S/N was not reached on AFGL 2136 because of forced closure of the telescope dome due to strong winds.

S140 IRS1 was observed on the same night from 08:53:38 to 09:05:19 and 09:45:41 to 09:57:22. BS8585 was taken as a standard star in between these acquisitions. BS8585, or α Lac, is an A1V star with V band magnitude of 3.77.

The data reduction was carried out using the Spextool package version 5.0.1 (Cushing et al. 2004). Wavelength calibration solutions were accurate to within 0.5 km s−1. The atmospheric seeing during the observations was good, therefore the binary in MonR2 IRS 3 was resolved, allowing for separate extraction of the two binary components. The spectra were then combined and the telluric correction was applied by dividing by the standard star, which removes the blaze function of the instrument and any fringing that was present.

The wavelength range of the L2 setting is in the middle of the H2O ice band. The final merged spectrum therefore had large scale structure in the baseline due to the ice feature. This was removed by fitting a baseline to the spectrum using a broad median filter and dividing the spectrum by this fitted baseline (Fig. 1). Some lines were affected by systematic error due to proximity to telluric residuals. In some cases a local continuum had to be used to achieve a flat continuum across the line by fitting a straight line across it. Many lines were lost in deep telluric lines, as the atmosphere is strong at these wavelengths.

Fig. 1 Full spectrum of the L2 setting of MonR2 1RS 3. The green line is a heavily smoothed median filter of the original spectrum used to define the baseline. The spectrum has been fully reduced, therefore the atmosphere and blaze function have been removed by division by the standard star. Regions of poor transmission have been masked out. The broad absorption feature at 3.041 µm is due to the presence of the H I (5–10) line in emission in HR1931, and the other broad absorption features are spurious features resulting from regions of very poor transmission.

3 Analysis

3.1 Emission lines

Emission lines are analysed following Barr et al. (2020), outlined in Sect. 5.3 of their paper. Consider a disk illuminating a (spherical) cloud of gas molecules at a distance d. The disk has a luminosity of 2L(v), where the factor 2 reflects the two surfaces of the disk. The continuum flux density seen at the Earth is then,

$$F_{\text{c}}=\frac{2L(v)}{4\pi D^2}$$(1)

where D is the distance to the Earth. If the gas cloud has a column density of N molecules, then the scattering intensity is,

$$I_{\text{s}}=N\frac{\text{d}\sigma }{\text{d}\Omega }\frac{L(v)}{4\pi d^2}$$(2)

with dσ/dΩ is the differential scattering cross section, which we will set equal to σ/4π. The observed flux density is then

$$F_l=\frac{\tau (v)}{4\pi }\frac{L(v)}{4\pi d^2}\frac{\pi r^2}{4\pi D^2}=\frac{\tau (v)}{4\pi }\frac{{\Omega }_{\text{cl}}}{2}{F}_{\text{c}}$$(3)

where r is the size of the scattering cloud or the size of the beam at the distance of the cloud, whichever is smaller, Ω_cl is the corresponding solid angle as seen from the continuum source, Ω_cl = πr²/4πd², and τ(v) is the optical depth at frequency v, τ(v) = Nσ(v). Here, we have made the assumption that d ≪ D. This solid angle represents a geometry factor. In the scenario where the scattering occurs in the upper layers of the disk, this geometry factor is ln[R_out/R_in] (Barr et al. 2020), resulting from an annulus in the disk.

Transitions were fit with a single Gaussian component with peak velocity, v_lsr, full width at half maximum (FWHM), and line depth as free parameters. Physical conditions were calculated using the rotation diagram method, where the column density in the upper level of each transition is given by:

$$N_{\text{u}}=8\pi \sqrt{2\pi }{\sigma }_{\text{v}}\frac{{\tau }_{\text{s}}}{A_{\text{ul}}{\lambda }^3}$$(4)

where A_ul is the Einstein A coefficient of the transition, λ is the wavelength, τ_s is the peak scattering optical depth, and σ_v is the standard deviation of the line in velocity units. We consider that the emitting gas fills the beam, but strictly speaking τ_s and N_u are a beam-averaged peak optical depth and column density, respectively. The FWHM is related to σ_v via the equation: $FWHM={\sigma }_{\text{v}}2\sqrt{2\ln (2)}$. For the H2 emission lines, a single Gaussian fit was not acceptable, therefore these lines are fit with a two Gaussian model with peak velocity, FWHM and line depth as free parameters. Average line profiles for low and high J lines of HCN are created by mathematically averaging normalised line profiles of each relevant line, interpolated onto the same velocity scale.

The above analysis depends on the assumption that emission lines form from resonant scattering of the infrared continuum, and also that the scattering gas is not in front of the continuum source and therefore the resonantly scattered lines are additive to the observed continuum flux. In the case that the emission lines arise from resonant scattering by molecules in a disk, the infrared continuum must be aligned with the plane of the sky, and the solid angle is contained in the geometry factor.

3.2 Absorption lines

3.2.1 Rotation diagram

Rotation diagrams for the absorption lines are created using the Boltzmann equation. The column density in the lower rovibrational level is calculated using the optically thin relation:

$$N_1=\frac{g_1}{g_{\text{u}}}\frac{8\pi }{A_{\text{ul}}{\lambda }_{\text{ul}}^3}{\displaystyle \int \frac{\tau (\upsilon )}{\varphi (\upsilon )}\text{d}\upsilon }$$(5)

N_l is the column density of the lower level, g_l and g_u are the statistical weights for the lower and upper level respectively, A_ul is the spontaneous emission coefficient for the transition, λ_ul is the wavelength of the line, τ(v) is the optical depth profile of the line in velocity space, and ϕ(v) is the normalised line profile.

3.2.2 Curve of growth for foreground absorption

We compare our rotation diagram analysis with a curve of growth analysis which considers a single absorbing slab of gas illuminated from behind by a background dust continuum source. The theoretical curve of growth in this way can be approximated by:

$$\frac{W_{\lambda }}{b{f}_{\text{c}}\lambda }~\frac{\sqrt{\pi }}{c}\frac{{\tau }_{\text{p}}}{1+{\tau }_{\text{p}}/2\sqrt{2}}$$(6)

for τ_p < 1.254 and:

$$\frac{W_{\lambda }}{b{f}_{\text{c}}\lambda }~\frac{2}{c}\sqrt{\ln [{\tau }_{\text{p}}/{\tau }_{\text{p}}/ln2]+\frac{{\gamma }^{\lambda }}{4b\sqrt{\pi }}({\tau }_{\text{p}}-1.254)}$$(7)

for τ_p > 1.254, taken from Tielens (2021) with the additional factor 1/f_c to take into account a covering factor. Here W_λ is the equivalent width, λ is the wavelength, τ_p is the peak optical depth, b is the Doppler width in km s−1 and γ is the damping factor. The expression for τ_p is given by:

$${\tau }_{\text{p}}=\frac{\sqrt{\pi }{e}^2{N}_1{f}_l\lambda }{m_{e}cb}$$(8)

where N_l and f_l are the column density and oscillator strength in the lower level, respectively, e is the electron charge, m_e is the electron mass and τ_p is the peak optical depth. In the empirical curve of growth, we introduce an additional parameter, f_c, the covering factor of the absorbing gas, which is divided into the equivalent width (e.g., the left hand side of Eqs. (6) and (7) become W_λ/bλf_c). This parameter describes the covering of the background illuminating source, where a smaller covering factor will result in a smaller observed equivalent width. The level populations are related by:

$$N_1=\frac{g_{1}N}{Q(T)}{\text{e}}^{-E_1/T}$$(9)

where Q(T) is the partition function, N is the total column density of the species, E_l is the energy of the lower level in Kelvin and T is a free parameter that will be of order of the temperature derived from the rotation diagram.

Adopting a column density, temperature, and Doppler parameter for the absorbing gas, the left-hand side of Eqs. (6) and (7), as well as τ_p, can be calculated for each line, resulting in an empirical curve of growth. Empirical curves of growth are created from the observed equivalent widths for a range of different values of T and N, and these are fit to the theoretical curve of growth. We also vary the Doppler parameter, as discussed in Sect. 5.4. The fit is quantified by a least square process where the best combination of T and N is chosen such that the reduced chi-squared is the lowest. The errors on these best fit parameters are taken from the 1σ contour.

3.2.3 Curve of growth for disk model

We also carry out a curve of growth analysis in the approximation of absorption lines being formed in a disk atmosphere, with a temperature decreasing with altitude in the disk, based on Barr et al. (2020). Curves of growth are constructed using the equation:

$$\frac{W_v}{2Y\Delta v}={\displaystyle {\int }_0^{\infty }\frac{{\eta }_{0}H(a,\upsilon )}{1+{\eta }_{0}H(a,\upsilon )}\text{d}\upsilon }$$(10)

where W is the equivalent width of the absorption lines in frequency space, ∆v is the Doppler width in units of Hz, and Y is the saturation point of lines and takes into account the temperature gradient in the atmosphere. H(a, v) is the Voigt profile as a function of the frequency shift, v (in velocity space and normalised to the Doppler parameter), and a the damping parameter (in Doppler units). The right-hand side of this equation can be calculated as a function of η₀, resulting in a theoretical curve of growth. We have used the simple approximations to this provided by Mihalas (1978). η₀ is the opacity at line centre set relative to the continuum opacity given by:

$${\eta }_0=\frac{{\kappa }_L(v=v_0)}{{\kappa }_{\text{c}}}=\alpha \frac{A_{ij}{\lambda }^3}{8\pi \sqrt{2\pi }{\sigma }_{\text{v}}}\frac{g_{\text{u}}}{g_1}\frac{X_1}{X_{\text{CO}}}$$(11)

where X_l is the abundance in the lower level. α is set by the requirement that the CO abundance, X_CO, is 10⁻⁴, and is essentially the dust opacity, σ_c, which results in σ_c of 7.7 × 10⁻²⁴ cm−2/g-of-gas, or α = 1.3 × 10¹⁹ for AFGL 2136 (Sect. 5.2 of Barr et al. 2020). The level populations are related by Eq. (9), however with the column density N replaced by the abundance with respect to CO X.

The definition of Y is given in the appendix of Barr et al. (2020) and is insensitively dependent on the temperature. The relevant temperature will be of the order of the derived from the rotation diagram. Curves of growth are considered for line formation from pure absorption (є = 1), pure scattering (є = 0) and a combination of absorption and scattering (є = 0.5), where e is a measure of the thermalisation of the line, such that the fraction of scattered and absorbed photons during line formation is (1 − є) and є, respectively (see appendix of Barr et al. 2020 for details).

4 Results

4.1 Structure of MonR2 1RS 3

MonR2 is a massive star forming region at a distance of 0.83 ± 0.05 kpc (Herbst & Racine 1976). In the star forming region are several protostars associated with different stages in the high mass star formation process. The dominant source at infrared wavelengths in the region is MonR2 IRS 3, which has a luminosity of 1.5 × 10⁴ L_⊙ (Henning et al. 1992). This object was itself resolved into three protostars in a deep bispectrum speckle interferometry study (Preibisch et al. 2002). MonR2 IRS 3A and MonR2 IRS 3B have masses of 12–15 M_⊙ and 8–12 M_⊙ respectively (Preibisch et al. 2002).

MonR2 IRS 3 has been imaged at 10.7 µm using the Keck telescope (Monnier et al. 2009). These observations reveal two sources, MonR2 IRS 3A and MonR2 IRS 3B. The mid-IR continuum of MonR2 IRS 3A originates predominantly from a non-spherical structure with a major and minor size of 144 ± 6 mas and 88 ± 5 mas respectively (~ 125 × 75 AU, with position angle (PA) of 41 ± 2°). Ten times fainter, extended emission in the north-east and south-west direction traces an IR reflection nebula, presumably associated with an outflow cavity. MonR2 IRS 3B lies to the north-east at a separation of 0.87″, and appears more like a point source, or a very small compact disk on the scale of <0.2″ (Monnier et al. 2009). Boley et al. (2013) find, in a very high spatial resolution interferometry study at 10.6 µm, that MonR2 IRS 3A is best described by non-spherical geometry with a major axis size of 173.0 ± 1.8 mas. The mid IR interferometry study had limited uv coverage and only a 1D Gaussian was fit to the data for MonR2 IRS 3B, resulting in a size of 37.6 ± 1.2 mas (Boley et al. 2013). MonR2 IRS 3A is clearly the dominant source of continuum with strong, extended nebulosity and K-band magnitude of 7.9 (Preibisch et al. 2002), consistent with the 10 µm relative fluxes (Monnier et al. 2009). MonR2 IRS 3B is also relatively bright in the K-band with a magnitude of 9.1, however it is less extended than MonR2 IRS 3A.

The presence of a disk in MonR2 IRS 3A is not well established with only hints of a disk (Yao et al. 1997; Preibisch et al. 2002; Monnier et al. 2009; Kwon et al. 2016; Fuente et al. 2021). Kwon et al. (2016) and Yao et al. (1997) find evidence for an optically thick disk in MonR2 IRS 3 from polarimetry studies at 2 µm where scattered light implied a non-spherical morphology. High degrees of linear polarisation are observed perpendicular to the elongation of the MonR2 IRS 3 nebula, hinting at a circumstellar disk oriented in the southeast-northwest direction. The spatial resolution of these studies did not distinguish MonR2 IRS 3A from MonR2 IRS 3B however. Therefore these results would support a disk in MonR2 IRS 3A, assuming that it is this source which dominates the scattered light in MonR2 IRS 3. MonR2 IRS 3A is itself extended in the northeast-southwest direction with a fan-shaped bipolar near-IR nebula, supportive of a geometrically thick disk or torus in the perpendicular direction (Preibisch et al. 2002). However, Boley et al. (2013) and Monnier et al. (2009) measure a position angles of 50° and 40° for the mid-IR continuum respectively, which would set this along the same orientation as the near-IR nebula, contrary to what would be expected from a disk perpendicular to the outflow direction. Thus the near-IR reflection nebula traces the inside of the outflow cone, and the 10 µm continuum emission traces hot dust in the same region.

HCN J = 3–2 sub-mm velocity maps, taken with the Atacama Large Millimetre Array (ALMA) reveal a blue-shifted, south-west extension, that coincides in position angle (60^°) and size (1650 AU) with the diffuse nebulosity around MonR2 IRS 3A seen in the mid-IR (Fuente et al. 2021). There is also red shifted emission lobe which lies closer to the central object (~ 1050 AU). As only one line has been measured and self-absorption effects are obvious, the temperature of the J = 3–2 emitting gas is not well determined (Fuente et al. 2021).

Emission from simple species such as SO2 and OCS, as well as more complex species such as CH3OH and CH3OCHO, was interpreted as a torus structure at ~175 K, encapsulating both

MonR2 IRS 3A and MonR2 IRS 3B. HCN was more spatially extended than all other species, ranging from −15 km s−1 to 22 km s−1, and potentially probing outflowing material, with more extended emission being filtered out (Fuente et al. 2021).

MonR2 IRS 3B is the source of a very collimated micro-jet with three knots seen in the K band along the north-east southwest direction (Preibisch et al. 2002), also indicative of of a disk in this source. Collimated jets and outflows are an indication that they are in an early stage of their evolution, where the jet or outflow becomes wider with time (Beuther & Shepherd 2005), implying that any disk associated with the jet of MonR2 IRS 3B would be in an early stage of its evolution and very compact, still undergoing a period of very active accretion. Pomohaci et al. (2017) detected H2 emission in the S(1) 2–1 line at 2.24 µm towards MonR2 IRS 3B, but not MonR2 IRS 3A, supporting the origin of these knots as shocked gas in a jet (Smith 1993; Eislöffel et al. 2000).

MonR2 IRS 3C is significantly weaker than, and spatially resolved from, MonR2 IRS 3A, thus we can be confident that the contribution of this third component to MonR2 IRS 3A is negligible. We do not spatially resolve MonR2 IRS 3C from MonR2 IRS 3B however as they are separated by only 0.25″. Also these two sources are comparable in terms of brightness (MonR2 IRS 3C is 0.5 mag fainter in K band), therefore it is possible that there is significant contamination from MonR2 IRS 3C, and as a result we do not draw any firm conclusions about MonR2 IRS 3B.

In Fig. 2, we illustrate the structure of MonR2 IRS 3. The extended nebulosity from MonR2 IRS 3A along the north-east south-west direction is shown, as well as the blue and red shifted components observed in CO in the mid-IR and HCN in the sub-mm. Yet to be observed circumstellar disks are indicated where they would be expected for both MonR2 IRS 3A and MonR2 IRS 3B, perpendicular to the outflow direction.

4.1.1 HCN emission lines

A total of 15 HCN lines from the v₁ band are detected towards MonR2 IRS 3A. All transitions are seen in emission, spanning an energy range of ~1000K. Also several low energy HCN emission lines are marginally detected towards MonR2 IRS 3B.

Average line profiles obtained through line stacking provide useful insights for carrying out the analysis. These are constructed for both high and low J lines, with the distinction being made at J = 4, and are shown in Fig. 3. We only use lines that are not affected by systematic error and are therefore more reliable. Fitting a single Gaussian to each line profile we find that the low J line profile is broader than the high J profile by 2.2 km s−1. Alsxo the peak velocity is slightly red-shifted, with 14.24 ± 0.04 km s−1 for the low J profile compared to 13.55 ± 0.02 km s−1 for the high J profile. The low J profile appears slightly non-Gaussian with a shoulder on the blue side. The systemic velocity of the cloud is 10.0 ± 1.0 km s−1 measured from submillimeter (sub-mm) observations (van der Tak et al. 2003). Therefore we fit an extra Gaussian to the low J profile and find that a weak contribution at 9 km s−1 fits well, likely reflecting cold gas from the cloud. In order to properly compare the high and low J lines, we therefore fit the low J lines (R1, R3, P3 and P4) with a two Gaussian model, fixing the parameters for the weak shoulder to those derived from the average line profile, and allowing the parameters of the main velocity component to be free. We then subtracted out this shoulder from the low J lines and re-calculated the stacked line profile, which is presented in Fig. 3. Fitting the re-calculated average line profile we find that the deconvolved line widths of the low and high J profiles come into agreement (7.7 ± 0.1 km s−1 and 7.9 ± 0.1 km s−1 respectively) however the difference in peak velocity remains, with 15.44 ± 0.04 km s−1 for the low J profile. This difference is a fraction 0.3 of a resolution element.

The contribution to the R(l) line at 00 km s−1 is very weak and only detected at 1σ. Therefore the error on this line, as well as the integrated strength of the main component, is very large. We calculate the error on this line by comparing it to the noise level at the same location in the spectrum. This is also the case for the P(4) line. We do not calculate the physical conditions of the component at 10 km s−1 due to the large uncertainties.

The rotation diagram of MonR2 IRS 3A is shown in Fig. 4 where the column densities of the low J lines are calculated from the line profiles after the contribution at 10 km s−1 was subtracted out. The R(1) line is excluded from the fit to the rotation diagram due to the large uncertainty in the removal of the weak shoulder. The data can be described by a single temperature component of 516 ± 71K and a column density of 1.6 ± 0.2 × 10¹⁶ cm−2.

We also tentatively present HCN in MonR2 IRS 3B with the R(l), R(3), P(3) and P(4) lines marginally detected. The deconvolved line width and v_lsr are 4.4 ± 1.2 km s−1 and 16.5 ± 0.5 km s−1 respectively. The average line profile is shown in Fig. 5 and shows only one velocity component. Correspondingly the rotation diagram indicates the presence of a single cold gas component at 48 ± 25 K and is presented in Fig. 4. The column density for this gas is 1.6 ± 0.8 × 10¹⁵ cm−2.

High spectral resolution observations of CO v = 1–0 taken with iSHELL are able to resolve the binary allowing for separate extraction of the two binary sources (Lee et al., in prep.). Line profiles from the contributions of each binary component are presented in Fig. 6. Line profiles reminiscent of P-Cygni profiles are observed in MonR2 IRS 3A, in both ¹²CO and ¹³CO, with absorption blue-shifted, and emission red-shifted, with respect to the systemic velocity. The red-shifted emission is centred at around 14 km s−1, in agreement with the HCN emission, and persists up until high J, being still present in the P(26) ¹²CO line, although significantly weaker than the P(14) ¹²CO line. CO and HCN line profiles are compared in Fig. 7.

Goto et al. (2003) observed ¹²CO v = 2–0 in absorption in MonR2 IRS 3 at a resolving power of 23 000 and found a two temperature component structure in the rotation diagram. In their study the binary components of MonR2 IRS 3 were not separated, but rather were analysed as one object. A temperature of 59 ± 4 K to 46 ± 3 K was found for the¹² CO cold component for Doppler parameters ranging between 2.5 and 4 km s⁻¹, respectively. The hot component of ¹²CO was approximately 250 K, not varying much with the chosen Doppler parameter, with a column density of 1.2 × 10¹⁹cm−2. The temperatures derived by Goto et al. (2003) for the main CO component, seen as blue-shifted absorption, are significantly lower than the HCN temperature we measure in emission for MonR2 IRS 3A. The hot CO temperature is however consistent with SO2 ($({225}_{-70}^{+50}\text{K)}$ K), also observed in absorption (Dungee et al. 2018).

The temperatures of the cold component of CO derived by Goto et al. (2003) are in agreement with the temperature of the HCN emission we derive for MonR2 IRS 3B. Since these temperatures are low, a natural conclusion would be to associate this gas with the envelope, however the velocity of HCN in MonR2 IRS 3B is red-shifted by 7 km s−1 with respect to the cloud velocity. The K band magnitude of MonR2 IRS 3C is only 0.5 mag fainter than MonR2 IRS 3B, and it is less extincted, therefore it could be possible that there is significant contribution from MonR2 IRS 3C included in the line profile for HCN in MonR2 IRS 3B, as these two sources are separated only by 0.25″. We do not draw conclusions on HCN in MonR2 IRS 3B due to the uncertainty in the results.

Fig. 2 Schematic outline of the potential environment of MonR2 IRS 3 where an expanding shell of gas is traced by different molecules at different wavelengths. The individual (yet to be observed) disks of MonR2 IRS 3A, MonR2 IRS 3B and MonR2 IRS 3C are given in yellow. The infrared reflection nebulae extending from MonR2 IRS 3 A are shown in black (Preibisch et al. 2002; Monnier et al. 2009). Scattering of continuum photons from MonR2 IRS 3A at 3 µm may produce the red-shifted HCN emission seen from the red side of an expanding shell, which is also evident in CO emission at 4.67 µm. The blue-side of the shell expanding towards the observer is only seen in CO absorption at infrared wavelengths, which is further away from the star and hence cooler. This side is also evident from HCN 3–2 emission at sub-mm wavelengths (Fuente et al. 2021). In this figure east is left and north is up.

Fig. 3 Line proflies created by averaging a selection of high quality low (bottom) and high-J (top) HCN lines toward MonR2 IRS 3A. J in this case is Jupper and the high energy lines are chosen as those with J ≥ 5, shown in blue, and the low energy lines are those with J ≤ 4, shown in orange. The vertical dashed line denotes the cloud velocity of 10 km s−1 (van der Tak et al. 2003). Panel a shows fits to the average line profiles given in black and individual contributions to the fit as blue dotted lines. Panel b shows fits to the average line profiles with the contribution at 9 km s−1 subtracted out.

4.1.2 H2 ro-vibrational lines

Two strong lines of H2 are detected in MonR21RS 3A, the O(4) and O(5) transitions, also in emission. These are clearly non-Gaussian and are fit with a model of two combined Gaussians; a narrow component superimposed on a broad underlying component (Fig. 8). For the narrow velocity component, the peak velocities of the 0(4) and 0(5) line are 3.4 ± 0.5 km s−1 and 5.4 ± 0.2 km s−1 respectively. The FWHM of the narrow component are 5.0 ± 0.3 km s−1 and 4.8 ± 0.2 km s−1 for the 0(4) and 0(5) lines respectively. The broad component has a FWHM and v_lsr of 40 ± 2.4 km s−1 and -5.0 ± 1.0 km s−1 respectively, for the O(5) line. The blue-shifted wing of the broad component in the O(4) line is affected by a telluric residual which removes the edge of the wing. Therefore we use the O(5) line as a reference since this line is in a clean part of the spectrum and not affected by systematic error, and fix the line width and velocity of the broad component when fitting the O(4) line to the values of the FWHM and v_lsr derived from the 0(5) line.

The same two lines of H2 are also detected in MonR2 IRS 3B and presented in Fig. 9. In this source, the lines are weaker, and three velocity components are present in the 0(5) line. The peak velocities of these components are −7.0 ± 5.0 km s−1, 3.5 ± 0.6 km s−1 and 13.5 ± 2.7 km s−1, with respective line widths of 20.4 ± 8.0 km s−1, 8.5 ± 1.9 km s⁻¹ and 14.1 ±4.4 km s−1. The most red-shifted velocity component is not apparent in the 0(4) line however, unless it has a significantly different velocity and has blended into the neighbouring component, creating the slight shoulder on the red side of this component. As in MonR2 IRS 3A, we fix the parameters of the velocity component in the 0(4) line that is affected by interference by an atmospheric residual to those of the 0(5) line. The peak velocity and line width of the central component is 1.8 ± 0.8 km s−1 and 11.8 ± 1.0 km s−1 respectively. The H2 lines likely trace the photon-dominated region (PDR) therefore we limit the discussion to HCN and do not analyse H2 further.

Fig. 4 Rotation diagrams of HCN lines in the v1 band detected in MonR2 IRS 3A (panel a) and MonR2 IRS 3B (panel b). The error bars displayed are the 1σ error from the Gaussian fitting, however the error on the column of the R(1) and P(4) lines is estimated from the noise level. Emission lines of which the error is underestimated from the fitting routine due to systematic error in the data reduction process are indicated in grey. These transitions are excluded from the fit used to derive the physical conditions. The column density of the low J lines in MonR21RS 3A is calculated using the integrated line strengths without the contribution from the cloud at 9 km s−1.

Fig. 5 Average line profile for lines detected in MonR2 IRS 3B. The vertical dashed line denotes the cloud velocity of 10 km s−1 (van der Tak et al. 2003).

Fig. 6 Line proflies from the 12CO and 13CO v = 1–0 band at 4.67 µm in the atmospheric M band (Lee et al., in prep.). The 12CO P(14), 13CO R(13) and 12CO P(26) lines are shown panels a, b and c respectively. The profiles in black and red originate from the 3A and 3B binary sources respectively. The dashed denotes the cloud velocity at 10 km s−1. The P Cygni profile is most obvious in the P(14) profile of the MonR2 IRS 3A source but actually present in all lines (Lee et al., in prep.).

Fig. 7 Normalised line profiles of CO, HCN and H2 from MonR2 IRS 3A. The orange and blue lines indicate the low and high J average profiles of HCN respectively. Finally the vertical dashed lines are at 10 km s−1 and 14 km s−1.

Fig. 8 Emission lines of H2 detected towards MonR2 IRS 3A. The global fit in indicated in red, whereas the individual Gaussian components that make up this are show as dotted blue lines. The vertical dashed line denotes the cloud velocity of 10 km s−1 (van der Tak et al. 2003). The feature at 35 km s−1 in the left panel is the HCN R(5) line.

Fig. 9 Emission lines of H2 detected towards MonR2 IRS 3B. The global fit in indicated in red, whereas the individual Gaussian components that make up this are show as dotted blue lines. The vertical dashed line denotes the cloud velocity of 10 km s−1.

4.2 AFGL 2136

4.2.1 HCN absorption

In contrast to MonR21RS 3, the HCN v₁ band detected in AFGL 2136 shows transitions exclusively in absorption. We detect 32 unblended lines of HCN, with many lines being blended with other lines from hot bands of HCN or C₂H2. We present the full spectrum of AFGL 2136 where absorption line are detected in Fig. A.6. Here, we restrict our analysis to the v₁ band, as our targeted S/N was not reached, and the weaker vibrationally excited bands are not strong enough for detailed analysis. The weighted average v_lsr and deconvolved line width of the v₁ band are 27.0 ± 0.5 km s−1 and 10.7 ± 1.3 km s−1 respectively.

A rotation diagram analysis results in a temperature of 540 ± 18 K and a column density of 1.4 ± 0.5 × 10¹⁷ cm⁻², and the rotation diagram is shown in Fig. 10. The four lowest energy lines deviate from a linear relation exhibiting a higher column than expected and are excluded from the fit. It is possible that the four lowest energy lines are slightly overestimated due to the presence of weak lines from a hot band of C₂H2 at the same frequency. This band is not detected in other lines, however, therefore it is very hard to remove. We cannot fit the four lowest energy lines to derive a temperature with a reasonable error due to the weak lines from hot bands, but it is also possible that these lines may be contaminated by cold foreground gas.

The measured line profiles for the HCN v\ band are consistent with those of the HCN v₂ = 2–0 band detected at 7 µm (Barr et al. 2020), and are presented in Fig. 11. Also the physical conditions derived from the rotation diagrams are in agreement between these two bands, with a temperature and column density of 592 ± 21 K and 1.8 ± 0.2 × 10¹⁷ cm⁻² respectively for the v₂ = 2–0 band. This contrasts to the v₂ band at 13 µm which, although having the same temperature and line profile, has a factor of 4 lower column density derived from the rotation diagram, despite all these bands probing the same lower level.

Therefore we carry out a curve of growth analysis to investigate the effects of high optical depth. Due to the uncertainty of the intrinsic line width, the curve of growth analysis is done for a range of different Doppler parameters, to explore how this parameter affects the derived physical conditions. We limit the Doppler parameter in such a way that it is less than the value consistent with the observed line width but greater than the thermal width, resulting in a range of values for 1.5 km s−1 < σ_v < 4 km s−1. Fits to the curves of growth are not good for σ_v < 3 km s−1 and there is not much difference between a choice of σ_v > 3 km s−1, therefore we elect to use the value that best fit the H20 absorption lines in AFGL 2136, which was a value of 3 km s−1 (Barr et al. 2022).

Initially we investigate a curve of growth analysis in the approximation of an absorbing slab of gas with covering factor, f_c, equal to 1 (see Sect. 3.2.2), for the v₁ band at 3 µm and the v₂v₂ = 2–0 band at 7 µm. Equivalent physical conditions are derived from the curve of growth compared to the rotation diagram. The curves of growth are presented in Fig. 12, where it can be seen that the lines have not yet transitioned onto the logarithmic part of the curve of growth.

From Fig. 13, we see that the absorption lines at 13 µm require the use of a covering factor in order to explain the observed equivalent widths. A covering factor of 0.3 fits the line profiles well. From a curve of growth with f_c = 0.3, we recover physical conditions that are in agreement with the curve of growth for the other bands. Without the use of a covering factor, the column density derived from the v₂ band is a factor of 4 lower, reflected in the column density difference derived from the rotation diagram analysis. Temperatures, column densities, and beam filling factors derived for the three bands are summarised in Table 2.

Therefore we also consider a curve of growth analysis following the disk model developed by Barr et al. (2020) to explore other potential locations of the absorbing gas in the proximity of the star. For the analysis we elect to choose the same Doppler parameter for the slab model. Overall the resulting temperature is not very sensitive to σ_v or є, however as σ_v increases, the best-fit abundance decreases, and the abundance is lower for є = 1 compared to є of 0 or 0.5. The curves with є = 1 give the best fits in terms of the reduced chi-squared, and these also give the best agreement for the temperature and abundance across the different wavelengths, thus we present these curves of growth in Fig. 14 and their results in Table 3. The curves of growth for the other values of є are given in the Appendix.

The column density difference with wavelength apparent in the rotation diagram analysis does not remain after carrying out the curve of growth for the disk model. The lines at 13 µm are beginning to transition onto the curve of growth for є = 1, and for є of 0 or 0.5, they are on the logarithmic part of the curve of growth (see Figs. A.2-A.5). The derived temperature of the three bands shows some variation however, as does the abundance, for a give value of є. Slight differences in the location of the continuum origin at the different wavelengths could produce some difference in the derived physical conditions between the two bands. Overall the three bands of HCN seem to trace the same physical conditions.

Fig. 10 Rotation diagrams of absorption lines in the HCN v1 band (panel a) and in the C2H2 v2 + (v4 + v5) band (panel b) in AFGL 2136.

Fig. 11 Average line profiles of the various bands of HCN (panel a) and C2H2 (panel b) detected in absorption towards AFGL 2136. The energy range over which the line profiles are averaged are given for each band. The dashed line represents the systemic velocity at 22 km s−1. The range of lower energy level E over which the average line profiles were calculated is shown for each band.

Fig. 12 Curves of growth for the v1 (3 µm), v2 v2 = 2–0 (7 µm) and v2 (13 µm) bands of HCN in AFGL 2136, calculated using the absorbing foreground slab model. Equivalent widths for the v2 and v2 v2 = 2–0 bands are taken from Barr et al. (2020). A covering factor of 1 is used for the 3 µm and 7 µm lines, and a covering factor of 0.3 is used for the 13 µm lines. Curves of growth are constructed for b = 4.5 km s−1 (σv = 3 km s−1 ).

Fig. 13 Absorption lines of varying opacity from the v2 band of HCN (top) and C2H2 v5 band (bottom) in AFGL 2136. The green solid line illustrates the expected line depth of the HCN and C2H2 lines at 13 µm, given the temperature and column density derived from the lines at 3 µm, and taking a covering factors of 1 (left) and 0.3 (right).

4.2.2 C₂H2 absorption

The v₂ + (v₄ + v₅) band of C₂H2 is detected towards AFGL 2136 at 3 µm. Nine lines are detected and all are in absorption. The peak velocity and deconvolved FWHM of this band are 26.4 ± 0.7 km s−1 and 7.5 ± 1.9 km s−1 respectively. From the rotation diagram (Fig. 10) we find a temperature and column density of 630 ± 235 K and 8.0 ± 3.1 × 10¹⁶ cm⁻² respectively. The line profiles are consistent with the C₂H2 (v₄ + v₅) band detected at 7 µm (Fig. 11), as is the temperature and column density. The temperature and line profiles of the v₂ + (v₄ + v₅) band are consistent with the 13 µm v₅ band, however the column density of the v₅ band is a factor of 7 lower than at 3 µm. Again, these three bands probe the same vibrational lower level.

For the curve of growth for the v₂ + (v₄ + v₅) band with the foreground slab model, we find temperatures and column densities that are consistent with the rotation diagram, and are derived using a covering factor equal to 1. These are presented in Fig. 15. Also, the lines are on the linear part of the curve of growth implying that the v₂ + (v₄ + v₅) band is optically thin and in LTE. These results are summarised in Table 2.

As for HCN, the C₂H2 lines at 13 µm require the use of a covering factor, as shown in Fig. 13. With a covering factor of 0.25, the temperature and column density of the v₅ band come into agreement with the other bands, and the line profiles agree with the LTE model. The lines have not yet transitioned onto the logarithmic part of the curve of growth. At 7 µm, no covering factor is needed and the column density for f_c = 1 is consistent with the 3 µm column, as is the line profile, although the error on the column is large. A covering factor of 0.5 would also give consistent results, therefore there is an error on f_c of the order a factor 2.

The three lowest energy lines of the (v₄ + v₅) band detected at 7 µm in Barr et al. (2020) are overestimated in the rotation diagram. This is also true for the curve of growth, and the fitting does not converge with these lines included. These lines are broader than the other lines of this band by 5 km s−1, therefore we leave them out of the curve of growth fitting, and attribute the larger line width to contamination by colder foreground gas along the line of sight. We also exclude the two lines that fall below the linear fit in the rotation diagram of the v₂ + (v₄ + v₅) band as these are affected by telluric residuals.

In the curve of growth for the disk model we again take σ_v = 3 km s−1. The best-fit curves for the v₂ + (v₄ + v₅), (v₄ + v₅) and v₅ bands are shown in Fig. 16, also for є = 1. For C₂H2, the curve of growth analysis only partially mitigates the column density trend observed in the rotation diagram, and we measure a temperature and abundance trend across the different bands with both increasing with decreasing wavelength (Table 4). This abundance and temperature trend from the curve of growth analysis is present for all values of є.

Fig. 14 Curves of growth for the v1 (3 µm), v2 v2 = 2–0 (7 µm) and v2 (13 µm) bands of HCN in AFGL 2136, shown in panels a, b and c respectively. Curves of growth are calculated using the disk model. Equivalent widths for the v2 and v2 v2 = 2–0 bands are taken from Barr et al. (2020). Curves of growth are constructed for a σv = 3 km s−1 and є = 1.

Fig. 15 Curves of growth for the v2 + (v4 + v5) (3 µm), (v4 + v5) (7 µm) and v5 (13 µm) bands of C2H2 in AFGL 2136, shown in panels a, b and c, respectively. Curves of growth are calculated using the absorbing foreground slab model. Equivalent widths for the (v4 + v5) and v5 bands are taken from Barr et al. (2020). A covering factor of 1 is used for the 3 µm and 7 µm lines, and 0.25 for the 13 µm lines. A Doppler parameter of 4.25 km s−1 is considered, consistent with the disk model (σv = 3 km s−1 ).

4.3 Orion BN, MonR2 IRS 2 and S140 IRS1

No molecular lines were detected in the spectra of Orion BN, MonR2 IRS 2 or S140 1RS 1, therefore we place upper limits on the HCN column density in these sources. Upper limits were calculated by comparing a model spectrum to the noise level at the expected location of the strongest lines for each source. Assuming a temperature of 500 K, we derive upper limits of <2 × 10¹⁵ cm⁻², <3 × 10¹⁶ cm⁻² and <6 × 10¹⁵ cm⁻² on the HCN column density for Orion BN, MonR2 IRS 2 and S140 IRS1 respectively.

Fig. 16 Curves of growth for the v2 + (v4 + v5) (3 µm), (v4 + v5) (7 µm) and v5 (13 µm) bands of C2H2 in AFGL 2136, shown in panels a, b and c respectively. Curves of growth are calculated for the disk model. Equivalent widths for the (v4 + v5) and v5 bands are taken from Barr et al. (2020). Curves of growth are constructed for a σv = 3 km s−1 and є = 1.

5 Discussion

5.1 HCN emission in MonR2 IRS 3A

5.1.1 Emission from expanding shell

The high temperature HCN emission from MonR2 IRS 3A is consistent in velocity with the P-Cygni emission feature of the CO profile in MonR2 IRS 3A (Fig. 6). Lee et al. (in prep.) do not calculate a temperature for the emission component of CO, however the presence of this emission feature in CO up to J = 26 implies that the CO gas is also hot. From the CO P-Cygni profiles, we infer that that there is an expanding shell of gas, where the blue-shifted side is seen in absorption against the continuum source, and the red-shifted side seen in emission traces the back side of the shell which is expanding away from the observer, illustrated in Fig. 2. The HCN emission would then be associated with the back side of this shell which is illuminated by the main continuum source in the region, MonR21RS 3A.

The HCN emission lines of MonR2 IRS 3B are consistent in terms of the line profiles of CO absorption in this source, and the HCN temperature is consistent with the cold component of CO observed in absorption by Goto et al. (2003) and Smith et al. (2020). The temperature is 48 K, and the column density is an order of magnitude lower than MonR2 IRS 3A, therefore we locate this as cold gas in the envelope of MonR2 IRS 3B. This low temperature gas creates ro-vibrational emission lines and therefore must be caused by resonant scattering from cold foreground gas.

We attribute the HCN high J ro-vibrational emission lines in MonR2 IRS 3A also to resonant scattering. We consider that this gas is heated by a shock responsible for the red emission lobe seen in the J = 3–2 map of HCN in the sub-mm. In principle, the emission lines could then represent thermal emission of HCN in the high temperature post shock gas. In a C-type shock, heating is over a long column density (10²¹ cm⁻²) as the ionised fluid frictions itself through the neutral fluid, reaching temperatures of around 2000–3000 K (Hollenbach et al. 2013). While a J shock is instantaneously heated and then cools down over a short column (10¹⁸ cm⁻²), heating due to H2 formation could keep the postshock gas at the order of 500 K over a column of 4× 10²²cm⁻² (Hollenbach et al. 2013). The critical density of the ro-vibrational transitions is 2× 10¹¹ cm⁻³ (Hernández Vera et al. 2017), however, some six orders of magnitude higher than for the pure rotational transitions. We consider such a high density unlikely and therefore analyse the emission lines assuming a resonant scattering process (Barr et al. 2020). However, to keep the HCN at 500 K via resonant scattering at 1000 AU from the star (the distance of the red-shifted HCN lobe apparent from the sub-mm) is challenging.

The temperature difference observed between the mid-IR CO absorption and HCN emission components in MonR2 IRS 3A likely reflects a temperature gradient that decreases as one moves away from the illuminating star. The temperature of the blue-shifted CO absorption is ~200 K, consistent with SO2 absorption (230 ± 15 K) which has comparable line profiles as CO (Dungee et al. 2018). The temperature of the absorption is comparable to the temperature of the sub-mm emission lines, ${175}_{-16}^{+15}$ K, which are described as coming from a circumbinary torus (Fuente et al. 2021). Thus the CO foreground absorption has expanded more than, and is further away from the star (~ 1650 AU) than the back side of the shell (~1050 AU), and hence cooler.

We place an upper limit on the blue-shifted absorption component of <2 × 10¹⁵ cm⁻² at a temperature derived from CO absorption, 200 K, an order of magnitude lower than the column of HCN seen in emission. Taking the ¹²CO column density from Goto et al. (2003), this corresponds to a HCN/CO abundance ratio of <1.7 × 10⁻⁴. We surmise that we do not see any HCN absorption at 200 K (foreground gas evident in the CO absorption lines) because the shell has expanded too fast and the gas is too cool to form HCN in high abundance (Bast et al. 2013; Agúndez et al. 2018). Also the emission may be extended - on the scale of the redshifted lobe (~1.3″) – so would appear brighter in the observed spectrum the more it fills the slit, whereas the absorption is only against the continuum source (~0.144″). The fact that CO is observed in absorption, but HCN is seen in emission, means that the CO gas is in front of, or mixed in with, the 4.67 µm continuum emission source, but that HCN is not in front of, or mixed in with, the 3 µm continuum source. Hence, if CO and HCN are both in a foreground cloud, this puts stringent constraints on the geometry as this cloud has to be very precisely positioned relative to the disk.

The shell would then correspond to a windblown bubble along the PA of the reflection nebula. The difference in peak velocity observed between the low and high J HCN lines suggests that these lines are probing the acceleration zone of the expanding shell or outflow. The low J lines probe deeper into the shell and therefore trace more red-shifted velocities compared to the high J lines which become optically thick closer to the source. Since the iSHELL L band spectrum was extracted at the continuum peak, if the HCN emitting gas is indeed at 1000 AU (1.3″) from the continuum peak, then the actual line flux density is conceivably much larger. This would imply a much higher line-to-continuum ratio and thus the derived column density is a lower limit.

For both binary components, we measure the scattering optical depth from the line-to-continuum ratio, which results in a column density (Barr et al. 2020; Sect. 5.3). This assumes that the 3 µm continuum flux density from the top of the disk (directly seen by us) is equal to the flux density from the bottom of the disk (seen by the back shell). The central source by far dominates the mid-IR continuum compared to the scattering nebulosity and therefore the line to continuum ratio will give a good estimate of the scattering optical depth. For a continuum normalised spectrum the optical depth is 0.01, therefore either 1% of the scattering photons are scattered into outline of sight, or the size of the scattering region is 1% that of the continuum. We cannot distinguish between the two of these. In this scenario where the HCN emission comes from a foreground cloud at distance, d, of 1000 AU, for a cloud size, r, of 500 AU, the column density of HCN from Eq. (3) would be 2.6 ± 0.3 × 10¹⁷ cm−2.

P-Cygni profiles are commonly observed towards massive YSOs in CO at mid-IR wavelengths (Mitchell et al. 1990, 1991). They have also been seen in CH4 towards the massive protostar NGC 7538 IRS 9 at 3 µm (Boogert et al. 2004), thus are a characteristic of these objects also in organic species. For CH4, the temperatures are similar to what we find for HCN, with 600 ± 100 K and 270 ± 100 K for the emission and absorption components respectively. Cooler CH4 could come from sublimation of ices (Boogert et al. 2004), however this is not to be expected for HCN which is not observed in ices (Theulé et al. 2013; Boogert et al. 2015).

The fact that P-Cygni profiles towards MYSOs seems common makes it less likely that it is fortuitous in MonR2 IRS 3A. Furthermore the lack of an observed disk along the direction perpendicular to the outflow makes the disk origin for HCN emission appear less likely, and the location of the HCN emission in the expanding nebula seems more attractive. However, the heating process of the HCN in the expanding shell scenario has challenges, whether by shocks or resonant scattering.

5.1.2 HCN emission from a circumstellar disk

Alternatively, the combination of blue-shifted absorption and red-shifted emission could be fortuitous in a complex environment of star formation where many different physical components could be present. In this case, the HCN emission lines in MonR2 IRS 3A may trace a disk photosphere instead of an expanding shell or outflow, as was inferred for AFGL 2591 (Barr et al. 2020).

In order to explain the results of a high resolution spectral survey in the mid-IR, we developed a model in which absorption lines probe a circumstellar disk around massive protostars AFGL 2591 and AFGL 2136 (Barr et al. 2020). Column density variations in HCN and C₂H2 of up to a factor of 10 were observed between ro-vibrational bands at 7 and 13 µm that trace the same lower level, with the bands at shorter wavelengths exhibiting a larger column. Ro-vibrational transitions of H2O are observed to saturate at non-zero flux, a direct prediction of the proposed model (Barr et al. 2022).

HCN emission lines were also detected towards AFGL 2591 at 3 µm, in contrast to the mid-IR spectrum where all species are observed in absorption (Barr et al. 2020). The temperature derived from the emission lines was 250 K compared to 700 K in the absorption lines, and the emission lines were a factor of 2 narrower than those seen in absorption. We interpreted the emission lines to originate from scattering of continuum photons seen under an angle from higher up layers of the disk photosphere, where temperatures are lower, compared to the absorption which probes deeper into the disk. Importantly, the presence of HCN and C₂H2 emission lines at 3 µm implied the absence of these species in the inner hot disk as otherwise we would see absorption lines (similar to the 7 µm and 13 µm HCN bands in this source).

Developing a similar scenario for Mon R2 IRS 3A, we would locate the HCN emitting gas high up in the disk photosphere where 3 µm photons from the inner disk are scattered into our beam. As the excitation temperature of the HCN emission lines is much higher than for AFGL 2591, this gas is closer to the central region. However, as for AFGL 2591, the region producing the 3 µm continuum photons cannot contain HCN. The fortuitous alignment of the scattering molecules discussed in Sect. 5.1.1 is therefore not an issue in the disk model, as HCN may be destroyed by OH radicals in the inner disk (Kress et al. 2008). In this case, the column density of HCN will be the same as that quoted in Sect. 4.1.1, namely 1.6 ± 0.2 × 10¹⁶ cm−2.

HCN emission may also potentially originate in the surface of an externally heated disk, analogous to low and intermediate mass stars (Carr & Najita 2008; Salyk et al. 2011). The luminosity of MonR2 IRS 3A is an order of magnitude lower than AFGL 2591, thus it is possible that the disk of MonR2 IRS 3A behaves more like those of lower mass protostars. Here, CO and SO2 mid-IR absorption (Goto et al. 2003; Dungee et al. 2018) would have to originate in the foreground, as while HCN would be present in the disk photosphere, CO and SO2 could not as this would produce emission lines. If MonR2 IRS 3A has an externally heated disk, HCN may be thermalised in the dense layers of the disk, and thus the emission be due to collisional excitation. However the details of such a model depend strongly on the dust and HCN temperature and optical depth distribution over the disk.

The mid-IR continuum of MonR2 IRS 3A is spatially extended along the PA of the outflow on scales of 0.2″. Therefore any disk that is present must be smaller than these scales (<170 AU) and is yet to be firmly detected. Multi wavelength observations of the HCN v₂ v₂ = 2–0 band at 7 µm and the v₂ band at 13 µm would allow detailed comparison of line profiles, peak velocities, temperatures and column densities. This might go far to elucidate the characteristics of the MonR2 IRS 3A system, and also the excitation mechanism of HCN emission. The spatial resolution of iSHELL is seeing limited to around 0.6″, therefore the location of the HCN emission lines relative to the continuum peak is something that could be mapped with iSHELL, which would help to understand the morphology of HCN in the infrared. Care should be taken then to avoid chopping on the array as that may lead to underestimated line flux density for extended emission.

5.2 Absorption lines in AFGL 2136

5.2.1 Foreground absorption

An optically thin LTE model grossly overestimates the observed HCN and C₂H2 lines at 13 µm, assuming a covering factor equivalent to the 3 µm and 7 µm lines. Therefore a much smaller covering factor is required for the 13 µm lines than for the other wavelengths. The results from the curve of growth analysis are in accordance with the rotation diagram analysis, which implies that the gas is optically thin. This is supported by the fact that all bands lie along the linear part of the curve of growth. The results for HCN compared to C₂H2 are very similar, which is to be expected due to the very comparable chemistry of these two species, with HCN/C₂H2 abundance ratios of ~1 commonly observed (Salyk et al. 2011; Agúndez et al. 2018).

The quality of the data does not allow us to determine the exact value for the covering factor at 3 µm and 7 µm, however we can constrain this parameter from the 13 µm band. For a covering factor f_c < 0.2, the fit to the 13 µm lines for both HCN and C₂H2 deteriorates, as the lines no longer lie along the linear part of curve of growth. This means that at 3 µm and 7 µm, the gas must at least almost fully cover the source, with f_c > 0.9 and f_c > 0.95 for HCN and C₂H2 respectively. This constraint on the covering factor of HCN and C₂H2 means that, in the case that these species are produced in a foreground cloud, these molecules must cover the source more than H2O, which has a covering factor of 0.4 at 6 µm (Barr et al. 2022).

In the comparisons between the bands at different wavelengths, we only consider the ortho states of C₂H2, since these are the only states detected at 3 µm and 7µm. Para states of C₂H2 were detected in the v₅ band at 13 µm however (Barr et al. 2020), therefore we can also carry out a curve of growth analysis for the para lines of C2H2 in this band. In doing so we derive a temperature of ${566}_{-87}^{+68}$ K, and a column density of ${4.1}_{-1.2}^{+1.1}\times {10}^{16}$ cm−2, taking f_c = 0.25. This results in an ortho-to-para ratio (OPR) of 1.6 ± 0.4, which is less than the expected value of 3 (Barr et al. 2020). Therefore, the curve of growth analysis and partial covering factors do not account for the low OPR observed in AFGL 2136. With the summation of the column density of both ortho and para states, we derive a total C2H2 column density of 1.1 × 10¹⁷ cm−2 for the v₅ band.

We conclude that the variation in the column density initially discerned from the rotation diagram analysis is due to dilution, not high optical depth effects in the 13 µm band. The coverage of the 3 µm and 7 µm bands may or may not be exactly the same, but it is very clear that the coverage of the 13 µm band is only ~0.3 of that at the shorter wavelengths. Thus the area of the 13 µm continuum must be ~3 times larger than the 7 µm continuum, but the HCN/C₂H2 containing foreground cloud must completely, or almost completely, “cover” the disk region responsible for the 3 µm and 7 µm continuum. We can estimate the size of the continuum emitting region, assuming a radial temperature gradient with an exponent in the range of −0.4 to −0.7 (Dullemond et al. 2007). Setting the temperature of the dust disk at 30 AU to 1200 K (the sublimation radius), we find that a temperature of ~500 K is reached at approximately 110 AU. This compares to the 10 µm size of AFGL 2136 of 400 AU (Monnier et al. 2009). High spatial resolution imaging of the dust continuum of AFGL 2136 at 3, 7 and 13 µm would help to assess the reality of the foreground absorption scenario.

5.2.2 Disk model

In contrast to AFGL 2591 and MonR2 IRS 3, HCN and C₂H2 in the L band of AFGL 2136 are observed to be in absorption. The line profiles are consistent with those detected from these species at 7 µm and 13 µm. Barr et al. (2020) proposed that the absorption lines from this source are coming from an internally heated disk, analogous to AFGL 2591. Abundance and column density variations between the 7 µm and 13 µm bands, which probe the same lower level, were interpreted as dilution of the longer wavelength lines due to a contribution of dust emission at the longer wavelength from a region in the disk with a very low HCN and C₂H2 abundance. These species require high temperatures to form in high abundances (Bast et al. 2013; Agúndez et al. 2018), therefore the high abundance material will be concentrated in the inner disk, and the outer colder parts of the disk do not contribute to the absorption lines.

In Sect. 4.3, we found that the column density variations of HCN could be explained by optical depth effects with the column density of the v₂ band at 13 µm being underestimated because it is optically thick. This is only the gas for absorption in the disk model. In the context of the circumstellar disk model, the fact that the HCN lines measure the same abundance implies that the HCN lines at 13 µm are not diluted. In this case, the high abundance HCN region in the disk must extend to the location of the 13 µm continuum origin. We note that an OPR of 3 can be achieved taking σ_v equal to 2 km s−1 and an є of 0.5, in the curve of growth analysis for the disk model. From Fig. A.5 it is clear that the C₂H2 lines are well onto the logarithmic part of the curve of growth. Thus, in this model, the low OPR of C₂H2 can be understood as a consequence of the o-C₂H2 lines becoming optically thick higher up in the disk photosphere compared to the p-C₂H2 lines.

While the abundance of the three bands of HCN are equal, the bands of C₂H2 at the different wavelengths show significant difference after optical depth effects have been accounted for with the curve of growth analysis: specifically, the abundance derived from the 3 µm lines is a factor 3 and a factor 6 higher than derived from the lines at 7 and 13 µm, respectively. This implies that the high abundance C₂H2 gas lies closer to the star compared to, and is less extended than, the high abundance HCN region. This scenario is supported by the higher temperature measured for C₂H2 at 3 µm compared to HCN (750 K compared to 450 K), however is surprising given the similarity in the chemistry of these species (Salyk et al. 2011; Agúndez et al. 2018). Taking the abundances from the 3 µm lines, we calculate a HCN/C₂H2 ratio of 0.68–0.76, depending on the value of є.

6 Conclusions

We have conducted a survey of five massive protostars from 2.95 µm to 3.25 µm to investigate how widespread transitions of organic species such as HCN and C₂H2 are in these objects, and how they relate to potential circumstellar disks. We detected transitions towards two out of the five objects, which are the sources in the sample known to exhibit hot core chemistry.

In MonR2 IRS 3, HCN transitions are detected in emission; confidently towards MonR2 IRS 3A and marginally towards MonR2 IRS 3B. A rotation diagram analysis shows that this HCN in MonR2 IRS 3A has a single temperature of 516 ± 71K and a column density of 1.6 ± 0.2 × 10¹⁶ cm−2. The line profiles of HCN are consistent with an emission shoulder of CO in this source and are red-shifted with respect to the cloud velocity. The CO emission feature is adjacent to blue-shifted absorption, indicative of a P-Cygni profile. The CO emission feature is visible up to high energy transitions, suggesting that the temperature is high, consistent with HCN.

The HCN emission of MonR2 IRS 3A could originate either from an expanding shell of gas, or from the atmosphere of a cir-cumstellar disk. The physical structure of HCN imaged in the sub-mm shows red-shifted emission at ~ 1000 AU along the orientation of the reflection nebula, most likely probing the outflow. If the high J ro-vibrational emission comes from the outflow cavity, resonant scattering of HCN due to 3 µm continuum photons is an unlikely mechanism to keep the gas at 500 K, and thermal heating due to shocks is also challenging due to the high critical densities of the ro-vibrational transitions. On the other hand, no disk has been confidently detected in MonR2 IRS 3A, and the ubiquity of P-Cygni profiles in high mass star forming regions makes it less likely that this profile is fortuitous. Therefore we cannot conclude for certain which of the proposed scenarios best fits out observations. Multi wavelength observations of HCN would help to clarify this, as absorption lines at 7 µm and 13 µm would be in strong favour for the presence of a disk, with emission lines at 3 µm probing the upper layers, analogous to AFGL 2591. Also mapping of the HCN emission with iSHELL would help to understand the location of the HCN emission lines relative to the continuum peak.

In AFGL 2136, HCN and C₂H2 transitions are detected, exclusively in absorption. Comparing the different bands observed at 3 µm, 7 µm and 13 µm which probe the same lower level, a rotation diagram analysis reveals variations in the derived column density. We therefore conducted two separate curve of growth analyses, assuming different models to explain the absorption lines; one a stellar atmosphere model to describe a circumstellar disk, and one an absorbing slab model to describe absorption due to a foreground cloud.

The curve of growth analysis for the absorbing slab model shows that, for this model, the lines of HCN and C₂H2 are optically thin across all wavelengths. In order to explain the observed equivalent widths, the lines at 13 µm require a partial covering factor of 0.3 and 0.25 for HCN and C₂H2 respectively. This compares to covering factors of >0.9 and >0.95 for HCN and C₂H2 respectively at both 3 µm and 7 µm, which have approximately equivalent covering factors. Thus, while the absorbing cloud must cover the 3 µm and 7 µm emitting area equally, 3/4 of the 13 µm continuum emission does not pass through this absorbing cloud. Based on the HCN and C₂H2 covering factors, these species must cover the source more than H2O, assuming a foreground origin. We derive an OPR less than 3, taking into consideration the partial covering factors.

For the disk model, the HCN abundance trend with wavelength previously observed can be attributed to optical depth effects as bands across all wavelengths measure the same abundance. For C₂H2, variations with abundance are observed after the curve of growth analysis, with higher abundances found at shorter wavelengths, consistent with a radial abundance gradient in a disk as proposed by Barr et al. (2020).

Therefore we cannot make a conclusion as to the location of the HCN emission in MonR2 IRS 3A. Multi-wavelength observations of this source would prove instrumental in further characterising the origin of hot HCN gas. In particular, the presence of much warmer gas in absorption at 7 and 13 µm – as for AFGL 2591 – would help solidify the presence of a warm disk, while dominance of absorption lines from gas with similar characteristics as CO would locate the absorption in the foreground.

References

Appendix A Additional material

Fig. A.1 Line profiles of the detected transitions of the HCN v1 band of MonR2 IRS 3A. Fits are shown in red, including a two Gaussian composite fit for the P(4), P(3), R(3) and R(l) lines with individual contributions show as blue.

Fig. A.2 The same as Figure 14 but with є = 0.

Fig. A.3 The same as Figure 14 but with є = 0.5.

Fig. A.4 The same as Figure 16 but with є = 0.

Fig. A.5 The same as Figure 16 but with є = 0.5.

Fig. A.6 L2 setting of AFGL 2136 in the range where absorption lines are present. In red is shown an LTE model of HCN at column density and temperature of 1 × 1017 cm−2 and 500 K respectively, and C2H2 at column density of 8 × 1016 cm−2 and 500 K respectively. The identifications of the HCN v1 band and C2H2 v2 + (v4 + v5) band are given as blue and green dashes respectively. The grey model is the atmospheric ATRAN model and regions of very bad transmission have been masked out of the science spectrum.

All Tables

All Figures

Fig. 1 Full spectrum of the L2 setting of MonR2 1RS 3. The green line is a heavily smoothed median filter of the original spectrum used to define the baseline. The spectrum has been fully reduced, therefore the atmosphere and blaze function have been removed by division by the standard star. Regions of poor transmission have been masked out. The broad absorption feature at 3.041 µm is due to the presence of the H I (5–10) line in emission in HR1931, and the other broad absorption features are spurious features resulting from regions of very poor transmission.

In the text

Fig. 2 Schematic outline of the potential environment of MonR2 IRS 3 where an expanding shell of gas is traced by different molecules at different wavelengths. The individual (yet to be observed) disks of MonR2 IRS 3A, MonR2 IRS 3B and MonR2 IRS 3C are given in yellow. The infrared reflection nebulae extending from MonR2 IRS 3 A are shown in black (Preibisch et al. 2002; Monnier et al. 2009). Scattering of continuum photons from MonR2 IRS 3A at 3 µm may produce the red-shifted HCN emission seen from the red side of an expanding shell, which is also evident in CO emission at 4.67 µm. The blue-side of the shell expanding towards the observer is only seen in CO absorption at infrared wavelengths, which is further away from the star and hence cooler. This side is also evident from HCN 3–2 emission at sub-mm wavelengths (Fuente et al. 2021). In this figure east is left and north is up.

In the text

Fig. 3 Line proflies created by averaging a selection of high quality low (bottom) and high-J (top) HCN lines toward MonR2 IRS 3A. J in this case is Jupper and the high energy lines are chosen as those with J ≥ 5, shown in blue, and the low energy lines are those with J ≤ 4, shown in orange. The vertical dashed line denotes the cloud velocity of 10 km s−1 (van der Tak et al. 2003). Panel a shows fits to the average line profiles given in black and individual contributions to the fit as blue dotted lines. Panel b shows fits to the average line profiles with the contribution at 9 km s−1 subtracted out.

In the text

Fig. 4 Rotation diagrams of HCN lines in the v1 band detected in MonR2 IRS 3A (panel a) and MonR2 IRS 3B (panel b). The error bars displayed are the 1σ error from the Gaussian fitting, however the error on the column of the R(1) and P(4) lines is estimated from the noise level. Emission lines of which the error is underestimated from the fitting routine due to systematic error in the data reduction process are indicated in grey. These transitions are excluded from the fit used to derive the physical conditions. The column density of the low J lines in MonR21RS 3A is calculated using the integrated line strengths without the contribution from the cloud at 9 km s−1.

In the text

	Fig. 5 Average line profile for lines detected in MonR2 IRS 3B. The vertical dashed line denotes the cloud velocity of 10 km s−1 (van der Tak et al. 2003).
In the text

Fig. 6 Line proflies from the 12CO and 13CO v = 1–0 band at 4.67 µm in the atmospheric M band (Lee et al., in prep.). The 12CO P(14), 13CO R(13) and 12CO P(26) lines are shown panels a, b and c respectively. The profiles in black and red originate from the 3A and 3B binary sources respectively. The dashed denotes the cloud velocity at 10 km s−1. The P Cygni profile is most obvious in the P(14) profile of the MonR2 IRS 3A source but actually present in all lines (Lee et al., in prep.).

In the text

	Fig. 7 Normalised line profiles of CO, HCN and H2 from MonR2 IRS 3A. The orange and blue lines indicate the low and high J average profiles of HCN respectively. Finally the vertical dashed lines are at 10 km s−1 and 14 km s−1.
In the text

	Fig. 8 Emission lines of H2 detected towards MonR2 IRS 3A. The global fit in indicated in red, whereas the individual Gaussian components that make up this are show as dotted blue lines. The vertical dashed line denotes the cloud velocity of 10 km s−1 (van der Tak et al. 2003). The feature at 35 km s−1 in the left panel is the HCN R(5) line.
In the text

	Fig. 9 Emission lines of H2 detected towards MonR2 IRS 3B. The global fit in indicated in red, whereas the individual Gaussian components that make up this are show as dotted blue lines. The vertical dashed line denotes the cloud velocity of 10 km s−1.
In the text

	Fig. 10 Rotation diagrams of absorption lines in the HCN v1 band (panel a) and in the C2H2 v2 + (v4 + v5) band (panel b) in AFGL 2136.
In the text

	Fig. 11 Average line profiles of the various bands of HCN (panel a) and C2H2 (panel b) detected in absorption towards AFGL 2136. The energy range over which the line profiles are averaged are given for each band. The dashed line represents the systemic velocity at 22 km s−1. The range of lower energy level E over which the average line profiles were calculated is shown for each band.
In the text

Fig. 12 Curves of growth for the v1 (3 µm), v2 v2 = 2–0 (7 µm) and v2 (13 µm) bands of HCN in AFGL 2136, calculated using the absorbing foreground slab model. Equivalent widths for the v2 and v2 v2 = 2–0 bands are taken from Barr et al. (2020). A covering factor of 1 is used for the 3 µm and 7 µm lines, and a covering factor of 0.3 is used for the 13 µm lines. Curves of growth are constructed for b = 4.5 km s−1 (σv = 3 km s−1 ).

In the text

	Fig. 13 Absorption lines of varying opacity from the v2 band of HCN (top) and C2H2 v5 band (bottom) in AFGL 2136. The green solid line illustrates the expected line depth of the HCN and C2H2 lines at 13 µm, given the temperature and column density derived from the lines at 3 µm, and taking a covering factors of 1 (left) and 0.3 (right).
In the text

	Fig. 14 Curves of growth for the v1 (3 µm), v2 v2 = 2–0 (7 µm) and v2 (13 µm) bands of HCN in AFGL 2136, shown in panels a, b and c respectively. Curves of growth are calculated using the disk model. Equivalent widths for the v2 and v2 v2 = 2–0 bands are taken from Barr et al. (2020). Curves of growth are constructed for a σv = 3 km s−1 and є = 1.
In the text

Fig. 15 Curves of growth for the v2 + (v4 + v5) (3 µm), (v4 + v5) (7 µm) and v5 (13 µm) bands of C2H2 in AFGL 2136, shown in panels a, b and c, respectively. Curves of growth are calculated using the absorbing foreground slab model. Equivalent widths for the (v4 + v5) and v5 bands are taken from Barr et al. (2020). A covering factor of 1 is used for the 3 µm and 7 µm lines, and 0.25 for the 13 µm lines. A Doppler parameter of 4.25 km s−1 is considered, consistent with the disk model (σv = 3 km s−1 ).

In the text

	Fig. 16 Curves of growth for the v2 + (v4 + v5) (3 µm), (v4 + v5) (7 µm) and v5 (13 µm) bands of C2H2 in AFGL 2136, shown in panels a, b and c respectively. Curves of growth are calculated for the disk model. Equivalent widths for the (v4 + v5) and v5 bands are taken from Barr et al. (2020). Curves of growth are constructed for a σv = 3 km s−1 and є = 1.
In the text

	Fig. A.1 Line profiles of the detected transitions of the HCN v1 band of MonR2 IRS 3A. Fits are shown in red, including a two Gaussian composite fit for the P(4), P(3), R(3) and R(l) lines with individual contributions show as blue.
In the text

	Fig. A.2 The same as Figure 14 but with є = 0.
In the text

	Fig. A.3 The same as Figure 14 but with є = 0.5.
In the text

	Fig. A.4 The same as Figure 16 but with є = 0.
In the text

	Fig. A.5 The same as Figure 16 but with є = 0.5.
In the text

Fig. A.6 L2 setting of AFGL 2136 in the range where absorption lines are present. In red is shown an LTE model of HCN at column density and temperature of 1 × 1017 cm−2 and 500 K respectively, and C2H2 at column density of 8 × 1016 cm−2 and 500 K respectively. The identifications of the HCN v1 band and C2H2 v2 + (v4 + v5) band are given as blue and green dashes respectively. The grey model is the atmospheric ATRAN model and regions of very bad transmission have been masked out of the science spectrum.

In the text