A Detailed View of the Circumstellar Environment and Disk of the Forming O-star AFGL 4176

The ALMA observations of AFGL 4176 were taken in 2014 August as a Cycle 1 Transfer project during ALMA Cycle 2. The observations are described in J15. Here, we mention several further properties of the observations. The cycle time of the phase calibrator was approximately 10 minutes, and the total time on-source across the observations was 1.37 hr. We number each of the four observed spectral windows (spw) by increasing frequency: spw0 to spw3. Table 1 provides further information about each of the observed ALMA spectral windows, including the frequency range covered, the imaged channel width (which is same as the observational spectral resolution), the imaged synthesized beam, and the image rms. The angular resolutions for each spw range from 026 to 030, corresponding to 1100–1300 au. Each of the four observed spectral windows were imaged using Briggs weighting and a robust parameter of 0.5, and the continuum was subtracted using imcontsub in CASA to avoid over-subtraction. The rms noise in each spectral window (given in Table 1) was measured after removal of the line emission, which was performed by σ-clipping each cube. Except for Figure 1, all figures in this paper show primary-beam-corrected images.

Zoom In Zoom Out Reset image size

Table 1.  Summary of Observed ALMA Spectral Windows

Spw Name	Frequency Range	Channel Width	Synthesized Beam	rms Noise
	(GHz)	(km s−1)		(mJy beam−1)
spw0	238.8376–239.3064	0.354	0.30 × 028 PA: 380	3.8
spw1	239.6035–241.4785	1.407	0.29 × 025 PA: −308	2.0
spw2	253.1055–254.9805	1.332	0.28 × 023 PA: −337	2.1
spw3	256.1146–256.5834	0.330	0.26 × 026 PA: −171°	3.4

Download table as:  ASCIITypeset image

When imaging H29α and C³⁴S, we instead used natural weighting to recover extended emission, producing beam parameters of 0.33 × 032 PA = 69.0° and 0.35×030 PA = −31.1°, respectively. The H29α image had an rms noise of 1.3 mJy beam⁻¹ in a channel width of 2.0 km s⁻¹ (larger than the native spectral resolution of 0.33 km s⁻¹), and the C³⁴S had an rms noise of 4 mJy beam⁻¹ in a channel width of 1.4 km s⁻¹. The ALMA and APEX C³⁴S images were combined using the CASA task feather.

We observed AFGL 4176 with the Australia Telescope Compact Array (ATCA) under project C2646 in the 15 mm band during 2012 on April 17, September 3, and December 21. The observations were designed to target transitions of NH₃ and hydrogen recombination lines. Table 2 presents a summary of the observations, including the start time, length of observation, range of baseline lengths, rms path length measured by the seeing monitor, largest angular scale (LAS), and primary beam for each of the three observing sessions. The pointing center of the observations was 13^h43^m019 −62°08'555 (J2000). Properties of the observed bands and lines on the three observing dates are given in Table 3. On the first two observing dates, only two continuum bands, with bandwidths of 2.112 GHz and comprising of 33 × 64 MHz channels, and two zoom bands, with bandwidths of 63.929 MHz and comprising of 2049 × 31.2 kHz channels, were available. This was expanded to two continuum bands and eight zoom bands on 2012 December 21; all of these zoom bands had a bandwidth of 63.929 GHz and 2049 channels, apart from the 23.709 and 24.509 GHz bands, which had bandwidths of 95.878 MHz and 3073 × 31.2 kHz channels.

Table 3.  Summary of Observed ATCA Bands

Observation	Central	Continuum	Line Rest	Imaged	Restoring	rms Noise	Detected?
Date	Frequency	or Line(s)	Frequency	Channel Width	Beam
(UTC)	(GHz)		(GHz)	(km s−1)		(mJy bm−1)
2012 Apr 17	20.434	Continuum	⋯	⋯	1.30 × 103 PA: 89	0.74	Y (Figure 16)
	23.712	Continuum	⋯	⋯	1.20 × 090 PA: 71	0.59	Y (Figure 16)
	20.466	H68α	20.46177	5	3.63 × 258 PA: 122a	1.85	Y (Figures 18, 19, 22)
	23.712	NH3(1, 1)	23.69450	0.4	12.64 × 1126 PA: 154b	6.51	Y (Figure 18)
	⋯	NH3(2, 2)	23.72263	0.4	12.67 × 1119 PA: 168b	6.34	Y
2012 Sep 2	24.139	Continuum	⋯	⋯	0.51 × 036 PA: 144	0.19	Y (Figure 17)
	24.533	Continuum	⋯	⋯	0.51 × 035 PA: 145	0.19	Y (Figure 17)
	24.139	NH3(4,4)	24.13942	0.8	0.54 × 038 PA: 148	1.95	Y (Figure 20)
	24.533	NH3(5,5)	24.53299	0.4	0.53 × 037 PA: 146	2.38	Y (Figure 20)
2012 Dec 21	20.905	Continuum	⋯	⋯	0.96 × 083 PA: −884	0.71	Y
	24.205	Continuum	⋯	⋯	0.87 × 077 PA: −887	0.59	Y
	20.457	H68α	20.46177	5	2.88 × 240 PA: −886a	2.74	Y
	21.385	H67α	21.38478	5	2.83 × 237 PA: −895a	4.42	Y
	23.405	H65α	23.40428	5	2.84 × 233 PA: −879a	3.93	Y
	23.709	NH3(1, 1)	23.69450	0.8	7.30 × 699 PA: 754c	14.25	N
	⋯	NH3(2, 2)	23.72263	0.8	7.25 × 699 PA: 747c	14.17	N
	23.885	NH3(3,3)	23.87013	0.8	7.20 × 709 PA: 327c	14.55	N
	24.141	NH3(4,4)	24.13942	0.8	0.91 × 084 PA: 862	7.08	N
	24.509	NH3(5,5)	24.53299	0.8	0.90 × 083 PA: 872	6.48	N
	⋯	H64α	24.50990	5	2.86 × 232 PA: −874a	3.18	Y
	25.069	NH3(6,6)	25.05603	0.8	0.88 × 081 PA: 854	6.30	N

Notes. Notes in the final column indicate the figures in which images of the observed band are shown.

^aA Gaussian taper of 2'' was applied to the data during imaging. ^bA Gaussian taper of 10'' was applied to the data during imaging. ^cA Gaussian taper of 5'' was applied to the data during imaging.

Download table as:  ASCIITypeset image

The setup calibrator, which was used to calibrate the delays, pointing, and an initial amplitude and phase calibration, was 0537-441 in all observations, and the primary flux calibrator was 1934-638. The bandpass calibrator was 1253-055, and the phase calibrator was 1352-63, which was observed in a cycle of 1 to 3 minutes with the target source. Calibration and imaging was carried out in the 20120419 ATNF version of MIRIAD. After flagging, the remaining continuum bandwidth was 1.8 GHz for each band in the 2012 April and September data sets, and 1.6 and 1.5 GHz for the 20.905 and 24.205 GHz bands observed on 2012 December. Due to poor weather on the final date of observation, we carried out self-calibration on the continuum data and applied these solutions to the line data, which improved the final images. The data were all imaged using a Briggs robust weighting parameter of 0.5. Unfortunately, none of the NH₃ lines observed in 2012 December were detected and, therefore, will not be discussed further.

We observed AFGL 4176 with the Atacama Pathfinder Experiment (APEX)⁸ 12 m antenna under program M0007_97 on 2016 April 23 and 26. The precipitable water vapor on these dates was 2.0–2.4 mm and 2.8–3.8 mm, respectively. The pointing center of the observations was 13^h43^m017 −62°08'513 (J2000). We used the Swedish Heterodyne Facility Instrument (SHeFI) APEX-1 receiver to observe C³⁴S(J = 5−4) at 241.016 GHz in the lower sideband, which covered frequencies between 240 and 244 GHz with a channel width of 76.3 kHz. The corresponding beam size at this frequency is 259. On-the-fly maps of $2^{\prime} \,\times$ 2' extent were taken in two perpendicular scan directions to reduce scanning artifacts. Pointing calibration was performed every 1 hr, the total observing time was 8.8 hr, and the total time on-source was 5.3 hr. Data reduction was performed within GILDAS/CLASS.⁹ The achieved rms in the final C³⁴S(J = 5−4) cube in the central part of the map is 66 mK in a channel width of 0.28 km s⁻¹ (three times the native resolution).

3.1.1. Millimeter Continuum Source Properties

In this section, we determine the observed properties of the 1.2 mm continuum sources detected with ALMA, as initially presented in Figure 1 of J15. In J15, only the observed properties of the brightest source mm1 were given; here, we present the properties of all of the continuum sources within 5'' of mm1. We measured the positions and flux densities of these sources using the astrodendro software.¹⁰ In this process, we used the non-primary-beam-corrected image and set the peak flux required for a detection to be 5σ, as the noise across this image is constant. The measured properties for these identified sources were then extracted from the primary-beam-corrected image. The flux limit above which flux densities are measured in the primary-beam-corrected image is 2σ in the non-primary-beam-corrected image, and the required flux difference between embedded structures to be assigned as a separate source is 1σ. Table 4 provides the properties of these sources, which include source name, peak position, peak flux S_peak, integrated flux S_int, mass M, peak column density ${N}_{{{\rm{H}}}_{2}}$, and the geometric mean diameter D, calculated from the area covered by each dendrogram structure. Figure 1 presents the non-primary-beam-corrected ALMA 1.21 mm continuum emission from AFGL 4176 as contours; the sources are labeled with the names given in Table 4, and the area of each source shown in grayscale. We detect 17 sources in the field. The peak and integrated flux densities measured for mm1 are slightly larger than that quoted in J15, which was determined using a Gaussian fit. This is mainly because the Gaussian source lies on a background of more diffuse emission, which was not included in the Gaussian fit but is included when the source fluxes are determined using astrodendro, which assigns groups of pixels to a source. The uncertainties on the flux densities given in Table 4 only give the random uncertainties, whereas an absolute flux calibration error of 10–20% is also expected.

Table 4.  Measured and Calculated Properties of the ALMA 1.2 mm Continuum Sources

Source	Peak Position	${S}_{\mathrm{peak}}$	${S}_{\mathrm{int}}$	Massa	Column Densitya	Diameter
Name	(J2000)	(mJy beam−1)	(mJy)	(M⊙)	(cm−2)	(au)
mm1	13:43:01.693–62:08:51.25	38.49 ± 0.08	61.22 ± 0.04	9.4b	8.3 × 1024b	3500
mm1 branchc	13:43:01.693–62:08:51.25	38.49 ± 0.08	86.21 ± 0.07	13.3b	8.3 × 1024b	6500
mm2	13:43:01.614–62:08:50.40	6.79 ± 0.08	11.87 ± 0.03	3.6	2.9 × 1024	2600
mm3d	13:43:01.886–62:08:55.35	2.07 ± 0.08	9.41 ± 0.05	2.8	8.7 × 1023	4200
mm4	13:43:01.886–62:08:52.65	2.05 ± 0.08	2.38 ± 0.03	0.7	8.7 × 1023	2700
mm5	13:43:01.764–62:08:53.80	1.54 ± 0.08	9.03 ± 0.05	2.7	6.5 × 1023	5100
mm6	13:43:01.714–62:08:50.20	1.30 ± 0.08	1.80 ± 0.02	0.5	5.5 × 1023	2100
mm7	13:43:02.100–62:08:52.65	1.21 ± 0.08	9.68 ± 0.06	2.9	5.1 × 1023	5400
mm7a	13:43:02.100–62:08:52.65	1.21 ± 0.08	2.76 ± 0.02	0.8	5.1 × 1023	2200
mm7b	13:43:02.042–62:08:53.40	1.06 ± 0.08	1.61 ± 0.02	0.5	4.5 × 1023	1800
mm8	13:43:01.993–62:08:52.20	0.94 ± 0.08	3.36 ± 0.04	1.0	4.0 × 1023	3500
mm9	13:43:01.536–62:08:50.10	0.71 ± 0.08	0.71 ± 0.08	0.2	3.0 × 1023	1300
mm10	13:43:01.365–62:08:50.65	0.66 ± 0.08	1.88 ± 0.03	0.6	2.8 × 1023	2600
mm11	13:43:01.179–62:08:50.40	0.65 ± 0.08	2.18 ± 0.03	0.7	2.7 × 1023	2900
mm12d	13:43:01.721–62:08:55.30	0.62 ± 0.08	1.13 ± 0.02	0.3	2.6 × 1023	2300
mm13d	13:43:01.557–62:08:55.15	0.60 ± 0.08	1.84 ± 0.03	0.6	2.5 × 1023	3100
mm14	13:43:01.514–62:08:51.00	0.60 ± 0.08	1.02 ± 0.02	0.3	2.5 × 1023	1900
mm15	13:43:01.308–62:08:50.40	0.57 ± 0.08	0.57 ± 0.02	0.2	2.4 × 1023	1400
mm16	13:43:01.265–62:08:48.50	0.53 ± 0.08	1.79 ± 0.03	0.5	2.2 × 1023	3100
mm17d	13:43:01.607–62:08:54.90	0.41 ± 0.08	0.57 ± 0.02	0.2	1.7 × 1023	1800

Notes.

^aThe masses and column densities are derived using opacities from Draine (2003a, 2003b). To obtain masses and column densities using opacities from OH94 for dust with thin ice mantles at densities of n = 10⁶ cm⁻³, the values should be divided by ∼4. ^bA temperature of 190 K was assumed, as determined in J15. ^cmm1 branch includes sources mm1 and mm2. ^dAs discussed in Section 3.1.1, these sources are likely to be dominated by free–free emission from ionized gas.

Download table as:  ASCIITypeset image

The mass M and peak column density ${N}_{{{\rm{H}}}_{2}}$ for each source were determined using

$$\begin{eqnarray}&&M=\displaystyle \frac{{{gS}}_{\mathrm{int}}{d}^{2}}{B(\nu ,T){\kappa }_{\nu }}\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&{N}_{{{\rm{H}}}_{2}}=\displaystyle \frac{{{gS}}_{\mathrm{peak}}}{2.8{m}_{{\rm{H}}}{\rm{\Omega }}B(\nu ,T){\kappa }_{\nu }},\end{eqnarray} \tag{ 2 }$$

where g is the gas-to-dust ratio, d is the distance, B(ν, T) is the blackbody function (a function of frequency ν and temperature T), κ_ν is the frequency-dependent opacity, 2.8 is the mean molecular weight, m_H is the mass of a hydrogen atom, and Ω is the ALMA continuum beam area. As in J15, we assume a distance of 4.2 kpc. We assume a gas-to-dust ratio measured for the solar neighborhood of 154 (Draine 2011) and a temperature of 100 K, except for mm1, where we use a temperature of 190 K, which was the average temperature of the disk derived from modeling CH₃CN line emission using CASSIS¹¹ (J15). To be consistent with J15, the dust opacity at 1.21 mm was assumed to be 0.24 cm² g⁻¹, which was determined from the Draine (2003a, 2003b) Milky Way dust model with R_V = 5.5. This is a model that fits the dust opacities for the ISM in star-forming regions well, and in the rest of this paper, we use this model when referring to ISM dust. In comparison, the Ossenkopf & Henning (1994, hereafter OH94) dust opacity at 1.21 mm for dust with thin ice mantles, subject to grain growth and with densities of n = 10⁶ cm⁻³, is 1 cm² g⁻¹, which, if used instead, would reduce the masses and column densities by a factor of ∼4. We could not use this opacity for the radiative-transfer dust-continuum modeling presented in J15, because the OH94 dust models do not have the associated scattering properties required as input to the radiative-transfer dust code Hyperion (Robitaille 2011). If dust grain growth has occurred in AFGL 4176, then we would expect higher opacities than typical ISM values. The higher densities associated with massive star formation would lead to fast grain growth, but massive stars have short formation timescales (the lifetime of the MYSO phase is ∼1–3 × 10⁵ yr; Kuiper et al. 2010, 2015, 2016; Davies et al. 2011; Kuiper & Hosokawa 2018). Boley et al. (2013) found evidence for an enhanced population of 1.5 μm grains in AFGL 4176, suggesting that grain growth has started, but this would not be enough to change the millimeter opacities significantly.

The sources mm3, mm12, mm13, and mm17 are likely dominated by free–free emission (see Section 4.1); thus, their calculated masses and column densities should be viewed with caution.

We determined the spectral index of mm1 between spw1 and spw2 (between 240.541 and 254.043 GHz) by fitting the continuum data in these spectral windows. We applied the phase self-calibration (i.e., no amplitude self-calibration) derived from the continuum data to the wide spectral windows spw1 and spw2 and selected the channels without line emission to estimate the continuum. Spectral windows spw0 and spw3 were not used, as they contained only a small number of line-free channels. We then found the integrated flux in each line-free channel within a 1'' diameter aperture around the mm1 peak and subsequently fitted these continuum channel fluxes by a power law, shown in Figure 2. The resulting spectral index is 3.36 ± 0.14, where the uncertainty is determined from the scatter in the data. To estimate other sources of uncertainty in our determination of the spectral index, we compared the spectral index derived during our calibration for our phase calibrator J1308-6707 (found to be −0.60 between spw1 and spw2) to that derived from the fluxes in the ALMA Calibrator Source Catalog¹² (−0.48 between 91.5 and 343.5 GHz with both frequencies observed on the same date), which provided an estimate of the uncertainty: −0.12. Thus, including the contribution from the scatter in the fit, the combined error in the spectral index is ∼0.2. To confirm our result, we also imaged the continuum emission in the same wide spectral windows using clean and nterms=2 and obtained a similar average spectral index of 3.35 over pixels within 015 of the continuum peak (where the estimated error in the spectral index in these pixels was <2).

Zoom In Zoom Out Reset image size

The measured spectral index of 3.4 ± 0.2 can be explained by dust emission, with a contribution of ν² from the Rayleigh–Jeans tail of the Planck function (expected at 1.21 mm for dust hotter than 30 K) and the remainder from the frequency dependence of the opacity ${\nu }^{\beta }$. In this case, we would infer β = 1.4 ± 0.2. ISM dust has values of $\beta \sim 1.6$ (Draine 2003a, 2003b), which is within 1σ of our measured spectral index. In the following section, we assess the contribution to mm1 from free–free emission at 1.21 mm.

3.1.2. Contribution of Free–free Emission to mm1

In Section 4.1, we find centimeter continuum emission associated with mm1, which cannot be explained by dust emission alone. Therefore, the millimeter continuum emission from mm1 may also contain free–free emission. The spectral index estimated between 1.21 mm and 1.23 cm is 1.56 ± 0.15 (Section 4.1), which might be explained by almost-optically thick free–free emission with the same spectral index. Assuming a hypercompact H ii region with a size of 2000 au is producing the centimeter continuum emission (Section 4.1), and the temperature of the H ii region is 10⁴ K, the electron density would have to be ≳5 × 10⁶ cm⁻³ to ensure that the turnover frequency was above 24.328 GHz (corresponding to 1.23 cm; Mezger & Henderson 1967). Although high, this would theoretically be possible, as the hydrogen molecular number density derived from molecular line fitting is >8 × 10⁷ cm⁻³ (J15), which is larger than this value. However, the spectral index at 1.21 mm is 3.4 ± 0.2, indicating that dust emission is present, and the morphology and extent of the emission at the two wavelengths is dissimilar, implying that the emission at each wavelength is not tracing the same material. Therefore, it is likely that both ionized gas and dust emission contribute to the continuum emission at both wavelengths, but that dust emission dominates at 1.21 mm and ionized gas at 1.23 cm.

We now derive the fraction of ALMA continuum emission due to dust. When combining the contribution from dust and ionized gas, the total flux density ${F}_{\nu }$ at a given frequency ν within the ALMA spws is given by

$$\begin{eqnarray}&&{F}_{\nu }={F}_{{\nu }_{0}}\left[f{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{{\alpha }_{d}}+(1-f){\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{{\alpha }_{\mathrm{ff}}}\right],\end{eqnarray} \tag{ 3 }$$

where ${F}_{{\nu }_{0}}$ is the total flux density at ${\nu }_{0}$, f is the fraction of ALMA continuum emission due to dust at ${\nu }_{0}$, and ${\alpha }_{d}$ and ${\alpha }_{\mathrm{ff}}$ are the spectral indices of the dust and ionized gas emission, respectively.

Starting from Equation (3), the spectral index can be expressed as

$$\begin{eqnarray}&&{\alpha }_{\mathrm{obs}}=\displaystyle \frac{d{\mathrm{log}}_{10}{F}_{\nu }}{d{\mathrm{log}}_{10}\nu }\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{{F}_{{\nu }_{0}}}{{F}_{\nu }}\left[f{\alpha }_{d}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{{\alpha }_{d}}+(1-f){\alpha }_{\mathrm{ff}}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{{\alpha }_{\mathrm{ff}}}\right].\end{eqnarray} \tag{ 5 }$$

As the ALMA spws cover a small range of frequencies, we assume that the observed spectral index ${\alpha }_{\mathrm{obs}}$ across these spws is constant. Therefore, by evaluating ${\alpha }_{\mathrm{obs}}$ at ${\nu }_{0}$, we find

$$\begin{eqnarray}&&{\alpha }_{\mathrm{obs}}=f{\alpha }_{d}+(1-f){\alpha }_{\mathrm{ff}},\end{eqnarray} \tag{ 6 }$$

and thus,

$$\begin{eqnarray}&&f=\displaystyle \frac{{\alpha }_{\mathrm{obs}}-{\alpha }_{\mathrm{ff}}}{{\alpha }_{d}-{\alpha }_{\mathrm{ff}}}.\end{eqnarray} \tag{ 7 }$$

Assuming an upper limit for the spectral index of dust emission of 3.6 (derived from the dust model used in J15), an observed (i.e., combined) spectral index ${\alpha }_{\mathrm{obs}}$ of 3.4 ± 0.2, and an upper limit to the ionized gas spectral index at 1.21 mm of 2, we find a lower limit to the dust emission fraction f of 87 ± 13%. Thus, dust emission dominates at 1.21 mm, even in the case of the largest possible contribution from free–free emission.

As we do not know the exact contribution to the emission from the ionized gas at 1.21 mm, and because the relative astrometry is not sufficiently accurate, and the lower resolution of the ATCA data precludes our knowledge of the morphology of the 1.23 cm emission, we do not attempt to remove the free–free contribution from the ALMA continuum data. However, we have determined it only contributes <13% and, thus, the majority of the emission arises from dust.

3.1.3. Comparison to Maser Positions

Figure 3 shows the positions of several maser species detected toward AFGL 4176, including four Class II methanol masers from Phillips et al. (1998), the position of ground-state OH masers at 1665 and 1667 MHz detected by Caswell (1998), and four water masers detected by Walsh et al. (2014). Although the spatially linear methanol maser group (red circles) appears to not be coincident with the continuum peak of mm1, the absolute positional uncertainty of the ALMA data (orange circle, ∼01 or ∼1/3× the beamwidth) would nearly allow the methanol maser group to be coincident with the peak of mm1. However, in Section 3.2.4, we show that the thermal CH₃OH emission is also offset to the east of the millimeter continuum peak of mm1, so that the maser emission is probably occurring in the areas of brightest thermal CH₃OH emission. The OH maser (cyan square) is coincident with an extension of 5σ emission to the southwest of mm1. There are similar but less-prominent extensions to the northwest and southeast. Given the orientation of the disk as determined by J15, OH masers could potentially be tracing the outflow cavity walls, which is similar to that seen in G35.20-0.74 (de Buizer 2006). Finally, although the positions of the water masers (blue triangles) are stated to be less accurate than the other masers shown, there is one maser slightly offset to the west of mm1 and another two masers that are similarly offset from mm2. This near-coincidence and similar offset from the two brightest millimeter continuum sources in the field may indicate that the water masers are in fact coincident with these two millimeter continuum sources.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.4. Comparison to Mid-IR Observations

In this section, we present the spectra and detected lines toward AFGL 4176, and we investigate the various categories of spectral and spatial morphology seen across the detected species.

3.2.1. Spectra and Detected Lines

3.2.2. Line Morphological Categories

By inspecting the channel maps for each line, we found that the lines could be grouped into several morphological types, which will be reviewed in detail in the following sections. These types are as follows: disk-tracing, blue-dominant, red-dominant, and outflow-tracing lines. There is also a fifth category labeled as "unknown", which covers lines that did not clearly fit into one of the above categories, often because the emission from these lines did not have a high signal-to-noise. Tables A1–A4 list the morphological group for each detected line; Table A5 summarizes the number of identified lines for each morphological type for each species.

In Figure 6, we compare the emission from the 10 brightest lines (eight for the red-dominant lines) that had a unique identification and had no other detected line within 10 MHz (∼12 km s⁻¹) for the three morphologies with compact emission: blue-dominant, disk-tracing, and red-dominant. Figure 6 visually shows the difference between these morphology groupings. The blue-dominant lines show a prominent bar of emission to the east of the peak of mm1 in the blueshifted and sometimes central channels. They also have fainter emission in the redshifted channel. The disk-tracing lines are more kinematically symmetric, showing a velocity gradient progressing from east to west. We will show in Section 3.2.3 that these are morphologically and kinematically similar to the upper K-ladder transitions of CH₃CN and, therefore, trace the disk. The red-dominant lines are more similar to the disk-tracing lines but are generally more extended in the redshifted channel, with most possessing an arc of emission reaching up to the northwest.

Zoom In Zoom Out Reset image size

To quantitatively determine the similarity and differences within and between the compact morphology groups, we calculated average reduced χ² values between the different lines. For each pair of lines, we calculated a χ² value using three channels at −2, 0, and +2 km s⁻¹ from the ${v}_{\mathrm{LSR}}$ (−52 km s⁻¹), considering pixels in a $2\times 2^{\prime\prime}$ field of view centered on the peak continuum position of mm1, as shown in Figure 6. To make the comparison uniform, we added noise to each of the three channels, so that all channels had a peak signal-to-noise of 15. We also only considered lines with a peak signal-to-noise that was above 15 in all three channels before adding noise. The reduced χ² was calculated by comparing these three channel maps separately to the channel maps of another line at the same velocities, and an average χ² for a pair of lines was then determined by combining the χ² over these three velocity channels.

With a reduced χ² for each pair of lines in hand, we computed an overall average reduced χ² (with associated error in the mean) between morphological types, including comparing each line morphological type with itself. The results presented in Table 6 show that the reduced χ² values when comparing the line morphologies with themselves are close to one, indicating similarity, as expected (if the images were identical, modulo noise, this would produce a reduced χ² of ∼1), but that the average reduced χ² is statistically significantly higher when the χ² is calculated from the comparison between line types, thus showing that the morphology groups are quantitatively different from one another.

Table 6.  Average Reduced χ² Values for Comparison between the Three Compact Morphological Types

	B	D	R
B	1.78 ± 0.05	3.33 ± 0.08	3.55 ± 0.07
D	3.33 ± 0.08	1.63 ± 0.05	2.14 ± 0.05
R	3.55 ± 0.07	2.14 ± 0.05	1.35 ± 0.05

Note. "B" stands for blue-dominant, "D" for disk-tracing, and "R" for red-dominant.

Download table as:  ASCIITypeset image

As the emission from the outflow-tracing lines is more extended and clearly different to the first three morphological types, we show these in Figure 7, which has a larger field of view of 4''. The outflow-tracing lines are similar in that they all have secondary extended emission near the position of mm2, not seen in the other lines, which is most prominent in the −2 km s⁻¹ channel. This emission extends away from mm1 in a direction perpendicular to the disk midplane in mm1.

Zoom In Zoom Out Reset image size

In terms of chemistry, the 10 O-bearing molecules we detect have a mixture of morphologies, but, for the most part, they have predominantly blue-dominant morphologies. In addition, we detect nine N-bearing molecules, all of which morphologically trace the disk (i.e., are disk-tracing lines), except H₂CCN, whose lines are too faint to determine a morphology. Lastly, there are five S-bearing molecules and one further molecule (CH₃CCH) detected, which have a mixture of morphologies.

Figure 8 displays the variation of the fitted peak velocity of the detected molecules as a function of the upper transitional level energy (left panel) and linewidth (middle panel), as well as their fitted upper transitional level energy as a function of linewidth (right panel). All molecules with a detection >20σ are shown. The different line morphological types (disk-tracing, blue-dominant, red-dominant, and outflow-tracing) are shown as points of different colors with error bars. Since the fitting algorithm used for the simultaneous fit with hundreds of free parameters did not provide errors for the fitted parameters, the errors were instead estimated by testing the fitting algorithm used with Monte Carlo simulations of synthetic data to determine the relationship between the error in the velocity and linewidth of each line and the signal-to-noise and linewidth of that line. The error bars are, for the most part, smaller than the data points.

Zoom In Zoom Out Reset image size

Each line type displays velocities or linewidths scattered around different values. A dashed line showing the mean velocity and linewidth for each line type is shown in each panel (ambiguously identified or possibly blended lines are shown as empty points). The mean velocities for the solid points for each line type are −51.64 ± 0.21, −52.44 ± 0.15, −52.88 ± 0.11, and −53.27 ± 0.10 km s⁻¹ for red-dominant, disk-tracing, outflow-tracing, and blue-dominant lines, respectively. The left and right panels of Figure 8 show no trend in the line velocities or linewidths with E_up, indicating that the velocities and linewidths are dominated by kinematics. The middle panel of Figure 8 shows that the line morphologies also cluster by linewidth, with the blue-dominant lines clustering around a mean linewidth of 4.27 ± 0.11 km s⁻¹, whereas the red-dominant lines cluster around a mean of 6.90 ± 0.34 km s⁻¹. The disk-tracing lines have the largest linewidths, with values clustering around a mean of 8.00 ± 0.22 km s⁻¹, which is not unexpected if they trace fast-rotating gas close to the (proto)star. The fact that the linewidths of the blue- and red-dominant lines are less than those of the disk-tracing lines could be explained if these lines are predominantly tracing the emission on the blue- and redshifted side of the circumstellar structure, respectively.

Figure 9 shows the peak positions within integrated zeroth-order moment maps that were made for all lines with detections >20σ. The zeroth-moment maps were made by integrating a spectral slab centered on −52 km s⁻¹ with a width of 20 km s⁻¹. The peaks of the zeroth-moment maps are distributed differently for each line type. In the case of the disk-tracing lines, the peaks line up close to the disk midplane, slightly shifted to the blueshifted, northeast side of the disk. Most of the blue-dominant lines peak far into the blueshifted side of the disk, ∼022 or ∼920 au from the continuum peak to the northeast along the disk midplane at PA ∼ 60° (J15). At this angular separation from the continuum peak along the disk midplane, these points also appear to show a bar-like distribution that has a PA perpendicular to that of the disk midplane. The peaks of the unblended and unambiguous lines in the zero moment maps for the red-dominant lines (Figure 9) lie in a tight cluster slightly southwest of the continuum peak (on the redshifted side of the disk).

Zoom In Zoom Out Reset image size

3.2.3. Disk-tracing Lines

In J15, we presented the observed CH₃CN J = 13–12 K ladder emission from AFGL 4176, which was well-modeled by a disk in Keplerian rotation. However, including CH₃CN, we find a total of 55 lines detected within the same observation that have a similar symmetric velocity gradient and, therefore, are also likely to be tracing the disk. These molecules include eight of the nine (excluding H₂CCN, which is too faint) N-bearing molecules listed in Table 5, and HCOOH (formic acid). Therefore, given that all of the molecules that include nitrogen that are bright enough to determine their morphology are disk-tracing, the presence of nitrogen in the molecules of these observations appears to indicate that it is a good disk tracer.

To quantify the similarity between the CH₃CN J = 13–12 K ladder lines, which we previously determined traced the disk via radiative-transfer modeling in J15, and the remaining lines that we refer to as disk-tracing, we calculated a similar combined χ² to that produced in Section 3.2.2. In this case, we compared the CH₃CN J = 13–12 K = 7 line that traces the inner disk to the non-CH₃CN lines (which included the vibrationally excited CH₃CN lines not analyzed in J15). The resulting average reduced χ² value is 1.31 ± 0.06, demonstrating the similarity of the non-CH₃CN lines (as well as vibrationally excited CH₃CN) to the CH₃CN lines known to trace the disk. If the CH₃CN K = 7 line is compared to all of the non-CH₃CN lines (excluding the vibrationally excited CH₃CN lines), the average reduced χ² value is 1.32 ± 0.10, confirming that excluding the vibrationally excited CH₃CN lines does not significantly change the result.

Comparing with the results of Bøgelund et al. (2019), we note that they find the velocity gradients seen in four of our eight disk-tracing lines: NH₂CHO or HC(O)NH₂, CH₃CN, C₂H₃CN, and C₂H₅CN, as well as in CH₃OCHO. We categorize several CH₃OCHO lines as disk-tracing and several others as blue-dominant; therefore, we find that some of these lines trace the disk.

Figure 10 presents zeroth- and first-moment maps of representative lines for each of the nine disk-tracing molecules listed above, as well as the brightest unblended line for CH₃CN v₈ = 1. The first-moment maps were created using a 5σ cut, except for the faint NH₂CN line, for which we used 3σ. The example lines shown in Figure 10 are shown in bold in Tables A1–A4.

Zoom In Zoom Out Reset image size

Similar to that found for CH₃CN, all of the lines shown in Figure 10 display a velocity gradient across the source with blueshifted emission in the northeast and redshifted emission in the southwest. In Figure 11, we present the PV diagrams of the lines shown in Figure 10, which also shows their similarity to the CH₃CN line, providing further evidence that this group of lines is indeed tracing the disk. For instance, each PV diagram in Figure 11 shows a velocity gradient from blue- to redshifted from east to west and resembles one of the simulated model PV diagrams for a Keplerian disk shown in J15. A few lines, such as HC₃N and NS, also clearly show a curve tending to higher velocities at smaller radii in the lowest contours in the top left and bottom right quadrants of the PV diagram, as expected in the case of Keplerian rotation.

Zoom In Zoom Out Reset image size

Both J15 and Bøgelund et al. (2019) found that more excited lines, which trace hotter gas close to the central source, have smaller extents and steeper velocity gradients. Therefore, it is likely that the spatially compact lines trace the disk while the extended lines (e.g., low K transitions of CH₃CN and HC₃N) may also trace part of the envelope. This is supported by Figure 12, where we plot the area of all pixels brighter than half the peak flux density in the zeroth-moment map for the disk-tracing lines against their linewidth, showing a trend of more compact emission for lines with larger linewidths.

Zoom In Zoom Out Reset image size

Due to the fact that it traces extended emission and is blended with a CH₃OH line on the redshifted side, HC₃N is the least clear-cut of the disk tracers. However, HC₃N evaporates off dust grains at a similar temperature to CH₃CN (Collings et al. 2004; Jaber Al-Edhari et al. 2017), the HC₃N J = 28–27 transition we detect has a similar critical density to that of CH₃CN J = 13–12 K = 3 (>1 × 10⁶ and 1.5 × 10⁶ cm⁻³ at 100 K, respectively), and the energies of the upper levels of the transition are similar (177.26 and 144.63 K). Therefore, assuming similar abundances, it is not surprising that the extent of the HC₃N emission shown in Figure 10 is similar to that of CH₃CN.

Several of the disk-tracing lines trace only the inner several hundred to thousand astronomical units of the source, namely NH₂CN, CH₃CN v₈ = 1, C₂H₃CN, and C₂H₅CN. These lines have reasonably high excitation temperatures (205.25, 607.71, 278.02, and 169.27 K); although, lower abundances of these molecules also probably play a role in their smaller extent.

Formamide or HC(O)NH₂ is particularly interesting, in that the lines show double-peaked line structures (see spw2 in Figures 5 and 13), which may indicate a lack of lower-velocity envelope emission from these lines, which would "fill in" the line profile close to the line center. In Section 3.2.7, we discuss how the majority of lines (including those for HC(O)NH₂) are likely to be optically thin, so optical-depth effects will not be important. However, missing flux issues due to interferometric filtering of emission at larger scales may complicate this interpretation. Nevertheless, the moment maps shown in Figure 10 show that, along with vibrationally exited CH₃CN and HCOOH, the HC(O)NH₂ emission displays a clear velocity gradient, is confined to the inner region of the disk, and is therefore a good disk tracer.

Zoom In Zoom Out Reset image size

3.2.4. Blue-dominant Lines

We find a total of 85 lines that are categorized as blue-dominant, and that molecules categorized as blue-dominant are predominantly oxygen-bearing. There are five molecules that only have lines with a blue-dominant morphology in Tables A1–A4, which are CH₂CO, CH₃OH, HCOCH₂OH, CH₃OCH₃, and H₂CS. There are four molecules that have lines that are designated mostly blue-dominant but with a few categorized as disk-tracing: CH₃OCHO, C₂H₅OH, CH₃COCH₃, and aGg'-(CH₂OH)₂. The only oxygen-bearing molecule that does not have lines that are categorized as blue-dominant is HCOOH. Conversely, H₂CS is the only line that is not oxygen-bearing that is categorized as blue-dominant.

Figure 14 provides an example of the morphology of a blue-dominant line, showing the CH₃OH 17(3,15) → 17(2,16) A⁺⁻ line with a rest frequency of 256.228714 GHz. As introduced in Section 3.2.2, this example methanol line, as well as the other blue-dominant lines, shows a blueshifted bar-like morphology that is perpendicular to the disk plane. In this example, this is especially noticeable in the −53.78 km s⁻¹ channel. At higher, redshifted velocities, the emission reaches around the peak continuum position toward the southwest, suggesting a possible ring or shell-type structure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Notably, the distribution of the brightest emission bears resemblance to the distribution of the four known Class II methanol masers in the region (Phillips et al. 1998), which are shown as crosses in Figure 14 and in the bottom left panel of Figure 9. The brightest maser component, component C in Phillips et al. (1998), shown as a slightly larger cross in Figure 14, has a velocity of −54.5 km s⁻¹, which is blueshifted compared to the systemic velocity (∼−52 km s⁻¹) and is slightly blueshifted compared to the mean velocity of the blue-dominant lines (−53.3 km s⁻¹). As seen in Figure 4(a) of Phillips et al. (1998), the four maser points are all blueshifted, with velocities between −56 and −52 km s⁻¹. Given the spatial and velocity coincidence of the masers and thermal emission, it is likely that these are tracing the same physical structure.

Due to the low inclination of the system (i ∼ 30° between the rotation axis of the disk and the line of sight; J15), it is unlikely that inclination and, thus, optical-depth effects would be causing the asymmetry of emission (see also Section 3.2.7). Instead, there may be something in the disk or circumstellar environment that is desorbing the oxygen-rich blue-dominant molecules from the dust grains at this position to increase their local abundance, such as a shock.

Previously, Bøgelund et al. (2019) found that some of the O-bearing species they detect toward AFGL 4176 peaked 02 offset from the continuum peak; however, they also state there are no large differences in the spatial distribution of N- and O-bearing species. In contrast, as discussed above, we find that the blue-dominant lines, which show an asymmetric morphology often peaking in the blueshifted side of the disk and are mostly O-bearing, have a very different morphology to that of the disk-tracing or N-bearing species. Further, the association of the blue-dominant line emission with the methanol masers also indicates that different physical processes are involved with the production of the emission from disk-tracing and blue-dominant lines, and, thus, N- and O-bearing lines. The difference in our findings could be attributed to the fact that Bøgelund et al. (2019) fit Gaussians to the zeroth-moment maps of each line and determine the peak position from these fits, whereas we inspected the channel maps for each line (e.g., Figure 14) to determine its morphological type. We also note that Bøgelund et al. (2019) find a different velocity gradient in the first-moment maps of CH₃OH, C₂H₅OH, CH₃OCH₃, and H₂CS compared to the molecules they found that trace the disk rotation. However, the channel maps of the blue-dominant lines (e.g., Figure 14) show that these molecules do show a similar velocity gradient to that of the disk-tracing lines, but the emission is often brighter and more extended on the blueshifted side of the disk.

3.2.5. Red-dominant Lines

3.2.6. Outflow-tracing Lines

3.2.7. Effect of Optical Depth on Line Morphologies

3.2.8. The Nature of mm2

In Figure 16, we present the ∼1'' resolution ATCA 22 GHz or 1.36 cm continuum emission observed in 2012 April. The emission shows a large H ii region extending ∼10'' to the south of mm1 and is very similar in morphology to that presented by Ellingsen et al. (2005) at 8.59 GHz with ATCA in the 6A array. The source mm1 itself is coincident with the northwest part of this extended emission, but there is no specific peak at this position. Ellingsen et al. (2005) measured an integrated flux density of 363 mJy at 8.59 GHz for the extended H ii region with the 6A array. The corresponding flux density of the emission we measure at 23.712 GHz with the 1.5B array, the image with the lowest noise, is 300 ± 30 mJy (assuming a 10% absolute flux calibration error). Given that these observations have similar beam sizes, this would indicate a spectral index α of $-0.19\pm$ 0.14, which is consistent with an optically thin extended H ii region.

Zoom In Zoom Out Reset image size

The higher-resolution 24.328 GHz or 1.23 cm continuum image taken with the 6A array in 2012 September is shown in Figure 17, along with the ALMA 1.21 mm continuum as orange contours. There is a bright bar of 1.23 cm emission in the south of the image, which is coincident with the millimeter sources mm3, mm12, mm13, and mm17. Therefore, these millimeter sources likely form part of the extended free–free emission from the aforementioned bar seen at longer wavelengths, and the masses and column densities given in Table 4 for these sources should be viewed with caution as they may be contaminated or wholly explained by ionized gas emission.

Zoom In Zoom Out Reset image size

There is a second bar of 1.23 cm continuum emission to the north, associated with mm4 and mm7 and pointing radially away from mm1, which forms the brightest emission in a faint arc seen at <3σ. The source mm1 lies on the edge of this arc, suggesting the centimeter emission could be tracing the ionized edges of a cavity formed by the feedback from mm1. An alternative interpretation would be that the radially positioned bar forms part of a jet from mm1, in which mm4 is a knot.

Finally, there is an unresolved 1.23 cm continuum source associated with mm1, which is shown in further detail in the inset panel of Figure 17. A Gaussian fit to the source gives peak and integrated flux densities of 1.25 ± 0.24 mJy beam⁻¹ and 1.33 ± 0.44 mJy, respectively. Given that the source is unresolved, its dimensions are <048 × 018 (PA ∼ 1445), corresponding to a physical size of $\lt 2000\,$ au $\times 760\,$ au. The integrated flux density from a Gaussian fit to the ALMA emission at 1.21 mm is 50 ± 4 mJy (J15). The spectral index between 1.23 cm and 1.21 mm is therefore 1.56 ± 0.15. As the smallest spectral index for the emission that can be produced by dust emission is α = 2 in the case of extreme grain growth (for which β would tend to 0), this confirms that the emission seen at 1.23 cm cannot be explained by dust emission alone and that some of the emission is due to ionized gas.

We now determine the minimum contribution of ionized gas to the 1.23 cm continuum, assuming that the spectral index of the dust does not change with wavelength. As found in Section 3.1.2, the spectral index observed at 1.21 mm is 3.4 ± 0.2. Since the free–free spectral index cannot be larger than 2, corresponding to optically thick emission, we know that the dust spectral index therefore has to be >3.4 ± 0.2. Assuming a dust spectral index of 3.4 ± 0.2 as a lower limit and that all of the 1.21 mm emission is from dust (as an upper limit), we obtain a contribution of $\lt {0.019}_{-0.007}^{+0.011}$ mJy from dust at 1.23 cm, meaning that the dust contributes $\lt {1.4}_{-0.5}^{+0.8}$% to the total flux at this wavelength. Therefore, the 3σ upper limit to the contribution from dust is 3.8%.

In Section 3.1.2, we determine that the ionized gas emission contributes <13% of the total flux at 1.21 mm, and above, we have shown that the ionized gas contributes >96% of the flux at 1.23 cm. Thus, the resulting ionized gas spectral index between 1.23 cm and 1.21 mm is <0.7.

Figure 18 presents the first- and zeroth-moment maps of the ATCA NH₃(1, 1) emission from the region surrounding AFGL4176 on parsec scales. To aid the detection of the extended emission present in the data, we imaged the NH₃(1, 1) line with a Gaussian taper of 10'', which was a compromise to ensure detection and reasonable resolution. As previously shown in Johnston et al. (2014), there is a clear velocity gradient across the NH₃(1, 1) emission, oriented from northeast to southwest, which is in the same sense as the velocity gradient seen in the molecular gas detected with ALMA on smaller scales (e.g., Figure 10). Thus, we may be observing the large-scale rotation of the envelope or toroid in which the AFGL 4176 disk lies. Alternatively, the toroid may instead be a filament with a velocity gradient (possibly due to infall) from which the AFGL 4176 mm1 disk has inherited its sense of rotation. However, only one low-surface brightness clump (AGAL 308.944+00.121; Urquhart et al. 2014) lies close to the ATLASGAL source associated with AFGL 4176 mm1 (AGAL 308.917+00.122). Therefore, the cloud is not clearly linked to another high-mass clump that would cause such a large-scale mass flow. At a flux level of 45 mJy km s⁻¹, the diameter of the zeroth-moment emission spans ∼50'' or ∼1 pc, giving a radius of ∼0.5 pc or 100,000 au.

Zoom In Zoom Out Reset image size

Figure 18 also shows the H68α emission in red contours (which we discuss in more detail in Section 5.2), with a white box delineating the area covered by Figure 22. The peak of the H68α emission is roughly coincident with the center of the bar of continuum emission seen at 1.23 cm (Figure 17). Therefore, the millimeter emission observed with ALMA, including AFGL 4176 mm1, also lies at the far blueshifted end of the toroid, where the velocities are >−51 km s⁻¹, which, in fact, agrees with the velocities of the molecular lines detected with ALMA (see, e.g., Figure 8). The peak position of NH₃(1, 1), −13^h43^m0147−62°09'025 (J2000), is offset 114 or 48,000 au from the position of AFGL 4176 mm1. This large offset could suggest that parts of the rotating envelope may be locally unstable to fragmentation and that fragmentation is occurring on the blueshifted side of the envelope. It could also suggest that, in addition to the compact H ii region and AFGL 4176 mm1, there may be other forming stars within the envelope. However, we detect no ALMA 1.21 mm continuum emission at the position of the NH₃(1, 1) peak (although, as this lies 75 from the pointing center, the sensitivity here is decreased by 20%).

We determined the temperature in the large-scale envelope by fitting the NH₃(1, 1) and NH₃(2, 2) emission. We fit Gaussians to the main component, and hyperfine components in the case of NH₃(1, 1), for each pixel in the NH₃(1, 1) and NH₃(2, 2) cubes and determined the rotation temperature ${T}_{\mathrm{rot}}$ via

$$\begin{eqnarray}&&{T}_{{\rm{rot}}}=\displaystyle \frac{-{T}_{0}}{{\rm{ln}}\left\{\tfrac{-0.283}{\tau }{\rm{ln}}[1-\tfrac{{S}_{{\rm{peak}}}(2,2)}{{S}_{{\rm{peak}}}(1,1)}(1-{e}^{-\tau })]\right\}}\end{eqnarray} \tag{ 8 }$$

(Ho et al. 1979), where ${T}_{0}\,\simeq$ 41.5 K, τ is the optical depth of the main line, and ${S}_{\mathrm{peak}}(1,1)$ and ${S}_{\mathrm{peak}}(2,2)$ are the fitted peak fluxes associated with each line. The kinetic temperature ${T}_{{\rm{k}}}$ was then found by iteration using the relationship given by Walmsley & Ungerechts (1983) and updated by Swift et al. (2005),

$$\begin{eqnarray}&&{T}_{{\rm{k}}}={T}_{\mathrm{rot}}\left(1+\displaystyle \frac{{T}_{{\rm{k}}}}{{T}_{0}}\mathrm{ln}\left[1+0.6\exp \left(\displaystyle \frac{-15.7}{{T}_{{\rm{k}}}}\right)\right]\right).\end{eqnarray} \tag{ 9 }$$

Only pixels where the peak flux density in the main line was $\gt 7\sigma$ in both transitions were fit. When fitting for the optical depth τ, we found that the fitted values were close to zero near the center of the toroid or clump, while they were higher at the edges. Given that this effect is likely due to the lower signal-to-noise at the clump edges, we instead chose to hold τ at zero during the fitting, assuming, therefore, that the NH₃ clump is optically thin.

We also accounted for the effects of depopulation of the lower two energy levels at higher temperatures by applying a correction factor determined by comparing Equation (9) with the more exact results of Maret et al. (2009). The resulting correction can be approximated by the relation (with both temperatures in kelvin):

$$\begin{eqnarray}&&{T}_{{\rm{k}},\mathrm{corr}}=1.09\,{T}_{{\rm{k}}}-1.\end{eqnarray} \tag{ 10 }$$

The resulting temperature map is shown in Figure 19, with the integrated H68α emission shown by black contours. The map shows a gradient in temperature, which increases from south to north, peaking at the northwestern edge of the clump. Although this temperature peak does not coincide with the position of mm1, it is interesting to note that it is reasonably nearby. Thus, this temperature gradient may be evidence of heating by mm1 and the compact H ii region in the north of the clump.

Zoom In Zoom Out Reset image size

We determined the total H₂ mass of the toroid or clump from the NH₃ emission using Equations (15) and (16) of Rosolowsky et al. (2008) in the optically thin limit and assuming an NH₃ abundance of ${10}^{-7.5}$ (Urquhart et al. 2015). The maps of the linewidth and the flux of the NH₃(1, 1) derived from the NH₃ line fitting described above were used to determine maps of the parameters ${\sigma }_{\nu }$ and the brightness temperature (which can replace ${T}_{{ex}}\tau$ in equation (16) of Rosolowsky et al. (2008) in the optically thin limit), respectively. Via this method, we found a map of the column density across the clump and find the total H₂ mass to be 9200 M_⊙.

The integrated ATLASGAL flux density at 870 μm is 17.32 Jy (Urquhart et al. 2014). Using Equation (1) and assuming a dust opacity at 870 μm of 0.43 cm² g⁻¹ (Draine 2003a, 2003b, with R_V = 5.5), a temperature of 35.7 K (the median temperature determined in the analysis above), and the remaining assumptions laid out in Section 3.1.1, we determine the mass of the associated dust and gas clump to be 5100 M_⊙. If we instead assume an opacity of 1.85 cm² g⁻¹ (e.g., Urquhart et al. 2013) derived from the opacities of the 10⁶ cm⁻³ thin ice mantle dust model in Ossenkopf & Henning (1994), we obtain a dust-plus-gas mass of 1200 M_⊙. The mass of the clump determined from the Draine opacities is within a factor of two of the mass determined from NH₃, confirming that there is a large reservoir of mass available for star formation around the AFGL 4176 mm1 disk and its surrounding millimeter sources.

Figure 20 shows the first-moment map of the NH₃(5,5) line observed with ATCA, overplotted with contours of the integrated emission of both NH₃(4,4) and NH₃(5,5), as well as 1.23 cm continuum. In contrast to NH₃(1, 1) and NH₃(2, 2), the emission from these more highly excited transitions lies close to mm1, marked by the large plus sign in Figure 20. The NH₃(5,5) emission is brighter than that of the NH₃(4,4) and peaks 02 (840 au) to the east of mm1, whereas the NH₃(4,4) emission peaks to the southeast. Correspondingly, the velocities shown in the NH₃(5,5) first-moment map are blueshifted (∼ −55 km s⁻¹) compared to most of the lines observed with ALMA (e.g., Figure 8). The upper energy levels of these two transitions lie at 201.092 and 295.942 K for NH₃(4,4) and NH₃(5,5), respectively. Although no trend in velocity is seen with the upper energy level for the lines detected with ALMA, it may be that the NH₃(5,5) emission is tracing a patch of gas hotter than ∼200 K within the blueshifted part of the circumstellar material that is not detectable in NH₃(4,4). It is unexpected that the hottest gas is not found at the peak position and velocity of the continuum source mm1, where heating from the young star would excite the gas. This observation may be indicating that some process, such as shocks, could be heating the gas in the blueshifted side of the disk. Thus, it appears that similar physical processes are exciting both the NH₃(5,5) and the blue-dominant lines, such as methanol, that are presented in Section 3.2.4.

Zoom In Zoom Out Reset image size

While other molecules have comparatively compact emission, there are three tracers detected with ALMA that show extended emission that reaches distances >2'' from mm1: CH₃CCH, H₂CS, and C³⁴S. Figure 21 presents the combined ALMA + APEX C³⁴S (5−4) channel maps, which display emission extended perpendicular to the position angle of the major axis of the disk (shown as a dashed line). There is an arc of emission to the southeast of mm1, seen most clearly in the −51.31 km s⁻¹ channel, coincident with the 1.21 mm continuum sources mm5, mm7, and mm8 (which, themselves, also follow the same arc; see Figure 1). In the same channel, there is also a U-shaped arc of emission pointing to the northwest of mm1, with mm1 positioned at the base of the U. Given the small linewidth of the C³⁴S (5−4) emission (5.8 km s⁻¹) compared to the velocities measured for the CO (3-2) outflow presented in J15 (13–16 km s⁻¹ from the v_LSR), C³⁴S is unlikely to trace the fast part of the outflow but instead a slower wide-angle wind or dense material in the outflow cavity walls. Interestingly, another example of where C³⁴S was found to be tracing a wide-angle outflow was the Class 0 protostar HH212 (Codella et al. 2014), a source that also has signatures of a Keplerian disk. Codella et al. (2014) find that the C³⁴S emission becomes more collimated at higher velocities, which can also be seen to the northwest of mm1 in Figure 21 at −52.71 km s⁻¹ and again more collimated at −54.12 km s⁻¹. They suggest that this demonstrates that the wide-angle outflow of HH212 has an onion-like structure, with higher-velocity material closer to the jet axis, which may also be what we are seeing hints of here for AFGL 4176.

Zoom In Zoom Out Reset image size

Assuming that the C³⁴S emission traces the walls or outer shell of a parabolic outflow described as $z={z}_{0}{(\varpi /{\varpi }_{0})}^{1.5}$, where z is the height from the x-y midplane, ϖ is the cylindrical radius, and z₀ and ${\varpi }_{0}$ are constants, we can find the half-opening angle of the structure. Taking the inclination of 30° found from our modeling of AFGL 4176 carried out in J15, we measured the size of the southeast outflow in the −51.31 and −49.90 km s⁻¹ channels by approximating the emission by an ellipse with a fixed aspect ratio (shown in Figure 21). This ellipse was measured to have semimajor and -minor axes of ∼2.2 × 19 or ∼9200 × 8000 au, respectively. The distance on the sky of mm1 to the center of the ellipse is ∼14. Assuming an inclination of 30°, the true distance between the center of the ellipse and mm1 is therefore z₀ ∼ 28 or 12,000 au at a cylindrical radius of ${\varpi }_{0}=9200$ au. From this, we could determine that the half-opening angle of the outflow or cavity wall at 1.5 × 10⁵ au is ∼19°. This is in reasonable agreement with the half-opening angle of 10° at 1.5 × 10⁵ au that was assumed for our previous modeling, as this provided a good fit to the SED.

Hydrogen radio recombination line (RRL) emission was observed and detected with both ATCA and ALMA. The zeroth- and first-moment map of the H68α line observed with ATCA in 2012 April is shown in Figure 22. We display this line as it was detected with the highest signal-to-noise of all of the observed ATCA RRLs. The H ii region traced by the H68α emission, like most of the 1.23 cm continuum, lies to the south of AFGL 4176 mm1. The first-moment map shows that there is an NNW-SSE velocity gradient (red to blueshifted, respectively) in the ionized gas on scales of 10,000 au. This does not agree with the velocity gradient seen in the NH₃ emission, which agrees more closely with that of the disk rotation seen in tracers like CH₃CN. It is also interesting to note that the central velocity of the H68α line (∼−44 km s⁻¹) is different to that of the molecular tracers, which lie around −52 to −53 km s⁻¹ (e.g., Figure 8). We fitted the mean H68α spectrum measured within a 1'' radius circular aperture centered on 13^h43^m01672–62°08'555 (J2000). The resulting flux, central velocity, and linewidth are 18.0 ± 1.1 mJy beam⁻¹, −44.4 ± 0.8 km s⁻¹, and 29.0 ± 2.0 km s⁻¹, respectively, with the stated uncertainties determined only from the scatter in the data. Interestingly, the emission becomes more redshifted closer to AFGL 4176 mm1. If the ionized gas is flowing from AFGL 4176 mm1, this would indicate that the redshifted side of the ionized flow has a higher redshifted velocity closer to the source compared to the systemic velocity of mm1 (v_LSR = 52–53 km s⁻¹). This velocity structure might be explained by a turbulent jet model (Arce et al. 2007) and the fact that the gas is ionized by radiative ionization or strong shocks in the outflow. Alternatively, the velocity gradient could be tracing dynamics driven by a separate source within this H ii region and/or by champagne flows due to varying density in the surrounding cloud. This seems the most likely explanation, as the bright bar of 1.23 cm continuum emission in the south of Figure 17 is coincident with the peak of the integrated H68α emission seen in Figure 22 and is centered on the velocity gradient within it, indicating that an independent source associated with this bar of continuum emission is probably producing the H68α emission.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We were only able to detect the remaining RRLs observed in 2012 December (H64α, H65α, H67α, and H68α) by finding their mean spectra over the area of sky covering the April H68α emission. We also fit these lines with Gaussians and determine their peak flux densities to be 13.3 ± 1.4, 12.9 ± 1.9, 14.6 ± 1.9, and 14.9 ± 1.5 mJy beam⁻¹, respectively, which are consistent to within 1σ. Similarly, the velocities of the lines are $-43.9\,\pm$ 1.3, $-47.4\,\pm$ 2.4, $-41.9\,\pm$ 1.9, and $-44.5\,\pm$ 1.2 km s⁻¹, and the linewidths are 25.5 ± 3.2, 33.7 ± 5.7, 30.5 ± 4.5, and 24.4 ± 2.9 km s⁻¹, which are consistent to within 2σ. As no trends can be deduced from these, the line properties are thus consistent with a single flux, velocity, and linewidth and, therefore, likely trace a similar range of densities within the H ii region.

We also detect H29α with ALMA. This H29α line emission is extended and faint, so that the resulting images are not useful. However, we measured a mean spectrum in a 5''(R.A.) × 35(decl.) ellipse centered on the same position as the circular aperture used for the ATCA RRLs. This resulted in a 9.3σ detection; the peak flux density, central velocity, and linewidth from a Gaussian fit are 0.58 ± 0.1 mJy beam⁻¹, −47.8 ± 2.1 km s⁻¹, and 25.7 ± 5.0 km s⁻¹, respectively. Although, due to large uncertainties, the linewidths of the H29α and H64, 65, 67, and 68α lines are consistent within errors, we calculated the probability distribution of electron density n_e for each pair of lines (i.e., H29α paired with one of the remaining lines) using Monte Carlo statistics and the method outlined in Keto et al. (2008). We then multiplied these distributions to find the combined probability distribution for n_e and found the 99.7th percentile provided an upper limit of ${n}_{e}\lt 6.6\,\times \,{10}^{5}$ cm⁻³, which, given the size of the emission, would be consistent with a compact H ii region.

We begin our comparison to other forming massive stars with rotating disk-like structures by presenting an updated version of the lower panel of Figure 14 shown by Beltrán & de Wit (2016) in our Figure 23. This figure plots the ratio of the freefall timescale (${t}_{\mathrm{ff}}$) to the rotational period (${t}_{\mathrm{rot}}$) against the gas mass of the structure (${M}_{\mathrm{gas}}$). The ${t}_{\mathrm{ff}}/{t}_{\mathrm{rot}}$ ratio compares the timescale for the gas to collapse under self-gravity in absence of rotation (${t}_{\mathrm{ff}}$) to the time for the structure to complete one rotation (${t}_{\mathrm{rot}}$). Thus, low values of ${t}_{\mathrm{ff}}/{t}_{\mathrm{rot}}$ indicate that the structure will only complete a small number of rotations before it is accreted on to the central star, so it is unlikely to be able to reach a stable configuration within this time.

Figure 23 has previously been used to separate the two main types of rotating structures observed around forming massive stars: (1) toroids, rotating flattened envelopes, or pseudodisks, which are massive (>100 M_⊙) large (∼10,000 au) rotating structures that are not in centrifugal equilibrium and lie at the bottom right of the diagram, and (2) disks in Keplerian rotation that exist over many rotations and lie toward the top left.

The properties of new or updated sources are given in Table 7 and are shown as blue points in Figure 23 (apart from AFGL 4176, which is shown as a purple star). The updated values have been taken from the papers listed in Table 7; however, some values remain as stated in Beltrán & de Wit (2016).

Table 7.  Source Properties for New and Updated Sources Displayed in Figure 23

Name	d	Lbol	Mgas	Rcont	Rvel	vrot	M⋆	tff/trot	References
	(kpc)	(L⊙)	(M⊙ )	(au)	(au)	(km s−1)	(M⊙)
AFGL 4176	4.2	1 × 105	0.91	870	1700	4.0	25	1.03	1
W3(H2O)E	2.0	2 × 104	8.8a	1150	1000	4.0	15	0.31	2
W3(H2O)W	2.0	2 × 104	6.5a	1150	1000	3.0	15	0.27	2
G328.2551−0.5321	2.5	1.3 × 104	0.15	250	800	9.0	13.5	3.90	3
G328.2551−0.5321 inner env.	2.5	1.3 × 104	4.7	1500	2000	4.5	13.5	0.55	3
G31.41+0.31	3.7	4.4 × 104	21a	1200	2400	7.5	38b	0.48	4,5
G23.01−0.41	4.6	4 × 104	1.5	2000	2000	3.5	20	0.76	6
IRAS 18162−2048 MM1	1.7	2 × 104	1.3	291	1700	3.4	14	0.73	7,8
G11.92−0.61 MM1	3.37	1 × 104	0.84	480	850	6.5	34	1.23	9
Orion Src I	0.415	1 × 105	0.1	50	80	13.0	15	2.18	10,11
IRAS 16547−4247	2.9	1 × 105	1.8	∼1000	1190	4.6c	20	0.70	12
G17.64 + 0.16	2.2	1 × 105	1.0	120	120	17.9	45	1.16	13
G339.88−1.26	2.1	4 × 104	2.8	530	530	6.0	12	0.49	14

Notes.

^aThese masses have not been recalculated and instead have been scaled from those given in the respective papers, as a varying temperature was used in their original calculation. ^bDue to the updated distance to the source, this mass is an upper limit. ^cThis velocity was scaled to ${R}_{\mathrm{vel}}$ using the model in reference 11.

References. 1: Johnston et al. (2015), 2: Ahmadi et al. (2018), 3: Csengeri et al. (2018), 4: Beltrán et al. (2018), 5: Beltrán et al. (2019), 6: Sanna et al. (2019), 7: Girart et al. (2017), 8: Girart et al. (2018), 9: Ilee et al. (2018), 10: Ginsburg et al. (2018), 11: Plambeck & Wright (2016), 12: Zapata et al. (2019), 13: Maud et al. (2019), 14: Zhang et al. (2019a).

Download table as:  ASCIITypeset image

As assumed in Beltrán & de Wit (2016), we calculate the freefall timescale ${t}_{\mathrm{ff}}$ as

$$\begin{eqnarray}&&{t}_{\mathrm{ff}}=\displaystyle \frac{\pi }{2}\sqrt{\displaystyle \frac{{R}_{\mathrm{vel}}^{3}}{2{{GM}}_{\mathrm{gas}}}},\end{eqnarray} \tag{ 11 }$$

and ${t}_{\mathrm{rot}}$ the rotation timescale as ${t}_{\mathrm{rot}}=2\pi {R}_{\mathrm{vel}}/{v}_{\mathrm{rot}}$. We use the radius measured from the line emission (${R}_{\mathrm{vel}}$) for the radius in all cases except for the case ${R}_{\mathrm{vel}}\lt {R}_{\mathrm{cont}}$ where we used the radius measured from the continuum (${R}_{\mathrm{cont}}$) for the calculation of ${t}_{\mathrm{ff}}$. This was to ensure that ${t}_{\mathrm{ff}}$ was evaluated at a radius that contained all of the mass of the structure. Apart from Src I, which has optically thick continuum emission, the masses determined from the dust-continuum emission ${M}_{\mathrm{gas}}$ for all of the sources in Table 7 and Figure 23 (including AFGL 4176 mm1) have been recalculated to be consistent with the set of assumptions stated in Beltrán & de Wit (2016), namely, that the gas-to-dust ratio is 100, the opacity at 1.4 mm is 1 cm² g⁻¹ (Ossenkopf & Henning 1994), and the dust opacity index β is 2. The temperature used to calculate the mass for each source was taken from the papers in Table 7.

Figure 23 also shows several dashed lines, which show the value of ${t}_{\mathrm{ff}}/{t}_{\mathrm{rot}}$ for structures in Keplerian rotation around stars with masses of ${M}_{\star }=0,10$ and 25 M_⊙, taking into account both the mass of the rotating structure and the star such that

$$\begin{eqnarray}&&{t}_{{\rm{ff}}}/{t}_{{\rm{rot}}}={\left[\displaystyle \frac{{M}_{{\rm{gas}}}+{M}_{\star }}{32{M}_{{\rm{gas}}}}\right]}^{1/2}.\end{eqnarray} \tag{ 12 }$$

Equation (11) assumes spherical symmetry, which, given the flattened nature of most of these structures, leads to an overestimate of ${t}_{\mathrm{ff}}/{t}_{\mathrm{rot}}$ (Beltrán et al. 2014). In addition, Equation (11) does not include the mass of the central star M_⋆, which could change the value of ${t}_{\mathrm{ff}}$ significantly in cases where ${M}_{\mathrm{gas}}\lesssim {M}_{\star }$. In fact, if M_⋆ is included in the calculation of ${t}_{\mathrm{ff}}$ in both of these equations, this has the effect of "flattening" Figure 23, so that the y-axis then becomes a measure of the "Keplerian-ness" of the structure; the dashed lines also become one horizontal line that shows the positions of structures in Keplerian rotation in the figure. We will further discuss such theoretically motivated variants of Figure 23 in Kee et al. (submitted); but for now, we show Figure 23 as presented in previous studies for a consistent comparison.

As can be seen in Figure 23, the new or updated sources (blue points) lie mostly above the previously known sources in the "disk" area of the diagram. Furthermore, several lie close to or above the upper dashed line for Keplerian structures around a 25 M_⊙ star, indicating that recent studies have started to push the discovery of Keplerian-like structures around massive forming stars to even higher stellar masses, into the regime of forming O-type stars. AFGL 4176 falls close to the dashed line representing Keplerian disks around a 25 M_⊙ star, in agreement with our previous modeling (J15).

Below, we expand on our brief discussion that was given in Section 1 on the variety of properties seen in disks around forming massive stars. While several sources typify the picture of a scaled-up version of a disk around a low-mass star (e.g., AFGL 4176, G11.92-0.61 MM1, IRAS 18162-2048 or GGD27 MM1, G023.01-00.41, and G339.88-1.26; J15, Ilee et al. 2016; Girart et al. 2018; Sanna et al. 2019; Zhang et al. 2019a), there are, in fact, a range of disk morphologies found toward high-mass stars. These can include multiple-disk systems (e.g., NGC 7538 IRS1, W3(H₂O), and W33A MM1; Beuther et al. 2017; Ahmadi et al. 2018; Izquierdo et al. 2018, in the latter case, fed by a large-scale accretion streamer), small disks (e.g., G328.2551-0.5321 and G17.64 + 0.16; Csengeri et al. 2018; Maud et al. 2019), and no clear evidence of disks (W51 and IRAS 18566+0408; Ginsburg et al. 2017; Silva et al. 2017).

Within the cluster surrounding the forming high-mass star GGD 27 MM1, Busquet et al. (2019) found a lack of disks close to MM1 (<0.02 pc) and that larger and more massive disks only existed at distances greater than 0.04 pc. They also found that the disks were, on average, smaller than more distributed environments such as Taurus. One explanation for this finding is that the rich cluster environment leads to truncation of these disks via more frequent interactions, which, along with photoevaporation, may also go toward the explanation of smaller and less massive disks in the Orion Nebula Cluster (Eisner et al. 2018). As massive stars are almost always found in clustered environments (e.g., de Wit et al. 2005), it is expected that interactions and accretion streams in clusters are important mechanisms for shaping disks (Bonnell et al. 2003; Bate 2018), whereas, despite their high EUV luminosity, disks around forming high-mass stars are found to be less affected by photoevaporation due to their high optical depth (Tanaka et al 2017; Kuiper & Hosokawa 2018).

As further observations of forming OB stars at high resolutions are carried out, the distribution of their disk properties, for instance, disk mass and radius, will act as an important discriminant between the importance of the cluster-scale environment in massive star formation. By definition, the cluster-scale environment does not play an important role for models of core accretion (McKee & Tan 2003; Krumholz & Bonnell 2009). If disk properties are found to be strongly affected by their cluster-scale environment, this suggests that a more holistic picture of star formation is required that takes into account not just the stellar core but the whole environment in which a massive star is formed.

High-mass star formation regions are known for their chemical richness, especially in complex organic molecules (COMs), which are defined by astronomers as molecules with $\geqslant 6$ carbon atoms (Herbst & van Dishoeck 2009), and this is indeed seen in our observations of AFGL 4176. Observations of regions such as Orion KL and Sgr B2 have often been at the forefront in the detection of these molecules, due, respectively, to their proximity and luminosity (e.g., Sutton et al. 1985; Cummins et al. 1986).

Segregation between different groups of species, such as that seen in our observations, has also been observed in Orion KL and many other massive star formation regions, most notably between N- and O-bearing species (Blake et al. 1987; Feng et al. 2015). Often, the N-bearing species are found to be more compact and closely associated with the hot cores, whereas the O-bearing species display more distributed emission (Blake et al. 1987; Öberg et al. 2013; Fayolle et al. 2015; Feng et al. 2015). There is a range of explanations put forward to account for this difference; a recent hypothesis by Suzuki et al. (2018) is that the N-bearing molecules are enhanced by reactions in the hot gas (>100 K) after the warm-up phase and evaporation of the icy dust mantles, as hydrogenation of these molecules into other species is not efficient at high temperatures. However, in contrast, Quénard et al. (2018) find that formamide, one of the disk-tracing and N-bearing species we detect, requires radical–radical grain-surface reactions as well as gas-phase reactions to reproduce the observed abundances (for further discussion, see Bøgelund et al. 2019). Thus, to-date, there is no definitive explanation of the chemical segregation of N-/O-bearing species.

As presented in Section 3.2.2, the observed species in AFGL 4176 (including 12 COMs) fall into four distinct morphological groups: disk-tracing, blue-dominant, red-dominant, and outflow-tracing lines. Except for HCOOH, all disk-tracing species contain nitrogen. Conversely, all molecules that we detect that contain nitrogen trace the disk (excluding H₂CCN, which is too faint). As discussed in Section 3.2.3 and can be seen from Figure 10, these species only trace the inner several hundred to thousand astronomical units of AFGL 4176, similar to the results of the previous studies mentioned above. In this way, AFGL 4176 follows the template of N-bearing species being associated with the hottest and densest regions of hot cores.

Given previous difficulties in unambiguously identifying disks around high-mass stars, only a few studies have been carried out to-date that begin to examine their chemistry. Isokoski et al. (2013) previously studied three high-mass YSOs with rotating disk-like structures and compared their chemistry to YSOs at the time thought to not contain disks, and they found no significant differences. However, several of the forming stars that were chosen as part of the comparison group that did not have disks have since been shown to contain (multiple-)disk systems (e.g., Maud et al. 2017; Ahmadi et al. 2018).

Specific high-mass YSOs with disks that have been previously studied in relation to chemistry include:

AFGL 2591 VLA3. Jiménez-Serra et al. (2012) found chemical segregation toward an MYSO known to harbor a disk, with the emission from three groups of molecules falling into three morphological types: (1) those peaking on AFGL 2591 VLA3 (H₂S and ¹³CS), (2) those that avoided the continuum position and were double-peaked (HC₃N, OCS, SO, and SO₂), and (3) CH₃OH, which presented a ring-like morphology. Their observations only covered one N-bearing molecule (HC₃N), which fell into the double-peaked morphological group. They explained the ring-like structure in CH₃OH emission by FUV photodissociation of this molecule within the inner regions of this MYSO. Whereas molecules such as H₂S and CS that peak at the central position can be reformed by gas-phase reactions, once destroyed methanol cannot be efficiently reformed in the gas phase (Garrod et al. 2008). This picture assumes spherical symmetry, but in reality, the disk will provide some degree of shielding from dissociation for CH₃OH. Therefore, it is unclear whether CH₃OH would be fully destroyed within several thousand astronomical units of the star.

NGC 6334 I(N) SMA 1b. Hunter et al. (2014) observed the emission from several different molecules toward this source, whose kinematics were found to be consistent with a rotating and infalling (sub-Keplerian) disk, including the N-bearing species HC₃N, CH₃CN, CH₃CH₂CN, and HNCO. However, in this source, many of the other observed species (such as CH₃OH) also trace the disk, showing, in this case, that there does not appear to be any strong chemical segregation of N-bearing species. Instead, the high-temperature transitions are seen to peak toward the source position, whereas the lower-temperature transitions had a double-peaked morphology, indicating a temperature gradient. However, one molecule that did not fit with this picture was HNCO, which had a compact morphology. Hunter et al. (2014) suggested this was the result of the destruction of larger parent molecules to form HNCO in the gas phase in the high-temperature inner regions of the MYSO. They also suggested that the disk morphology could provide shielding from dissociation for molecules such as HC₃N, explaining its emission at high velocities close to the source.

G11.92-0.61 MM1. Ilee et al. (2016, 2018) observed a range of species toward this MYSO, finding a Keplerian disk with an enclosed dynamical mass of 40 ± 5 M_⊙. Species that traced the disk include CH₃CN, CH₃CH₂CN, HNCO, DCN, OCS, H₂CO, CH₃OH, and CH₃OCHO. However, although they displayed the same velocity gradient as other tracers, the emission from the observed methanol lines was offset to the southeast compared to the dust continuum. Many of the observed molecules displayed double-peaked and/or asymmetric emission, except for OCS, which was coincident with the central source. In their modeling of the CH₃CN emission, they find evidence for two temperature components, which they suggest may originate from the hot (>150 K) and dense gas near the midplane, where CH₃CN has thermally desorbed from ice mantles, and from the atmosphere of the disk where gas-phase formation of CH₃CN is more important.

G345.4938 + 01.4677 / IRAS 16562-3959. Guzmán et al. (2018) analyzed emission from 22 different species toward the MYSO G345.4938 + 01.4677, which has previously been found to contain a compact rotating core seen in SO and SO₂ (Guzmán et al. 2014). They categorized their observed molecules into two groups. The first group (the "Shock group") is comprised of species whose emission was similar to SiO, and the second group (the "Continuum group"), containing a broad range of species including COMs such as CH₃C₃N, is associated with the 3 mm dust-continuum emission. However, they note that except for SO and SO₂, no other molecules display the previously seen velocity gradient that indicated that a disk may be embedded in the rotating core.

G35.20-0.74 N. Allen et al. (2017) found that N-bearing species, specifically cyanides, were found to be more abundant in part B3 of the source G35.20-0.74 N B, which has previously been found to be a Keplerian disk by Sánchez-Monge et al. (2014). They suggested that this chemical segregation could be from fragments or different sources within a disk or rotating torus, where hot gas-phase chemistry is more active in source B3, producing the larger abundance of N-bearing species.

GGD27 MM1 (HH80-81). Girart et al. (2017) found that sulfurated molecules such as SO₂ and SO, as well as H₂CO traced a rotating disk. They attributed the production of the sulfurated molecules to shocks occurring at the radius of the centrifugal barrier, also seen toward low-mass YSOs (Sakai et al. 2014; Oya et al. 2016). Another explanation they put forward was UV photodissociation of water in the disk, which allows the reactants to form molecules such as SO₂. In contrast, they found that methanol departs from the emission expected for a rotating disk and that it instead traces the outflow cavity walls. Detected isotopologues of HCN and HC₃N also did not show a clear velocity gradient; however, this may have been due to the fact these detections had a low signal-to-noise.

Orion SrcI. Ginsburg et al. (2018, 2019b) investigated the dust continuum and molecular lines detected at three frequencies between 0.87 and 3.0 mm toward the disk and outflow in this source, which, at a distance of 415 pc, constitutes the nearest example of a massive star with a disk. They detected SiO, water, and a host of unidentified lines, which they went on to identify as isotopologues of NaCl, KCl, and possibly AlO. By fitting the position–velocity diagram of the H₂O 5_5,0 − 6_4,3 line, which they find is tracing the upper envelope of the disk and the lower section of the outflow, they determine a central mass of 15 ± 2 M_⊙. They find that the salt lines also trace the base of the outflow or an upper layer of the disk but one that lies closer to the disk midplane than water and non-vibrationally excited transitions of SiO. This interpretation is complicated by the fact SrcI is nearly edge-on, and the dust emission is optically thick (Plambeck & Wright 2016), masking any line emission close to the midplane. Given the rarity of salt lines in the ISM, Ginsburg et al. (2019b) suggest that they may be a unique tracer of disks around massive stars, yet, at least in SrcI, they do not trace the midplane but instead an upper layer in the disk or the base of the outflow.

G17.64 + 0.16. This MYSO was one of six sources studied by Cesaroni et al. (2017), who found it was the most evolved in their sample but also the most chemically rich. However, their analysis found clear evidence for a disk. Later, G17 was revisited by Maud et al. (2018) and at higher resolution in Maud et al. (2019). Maud et al. (2018) found a small (∼200 au radius) rotating disk in SiO, with a similar position angle to the continuum emission. Comparable to the results for SrcI (Ginsburg et al. 2018), the presence of SiO in the disk around this source indicates the presence of possibly ionized and turbulent hot shocked gas. Indeed, the detection of compact H30α emission toward the source supports this. Both CH₃CN and CH₃OH are not found to trace the disk but instead the cavity working surfaces at the point the wide-angle wind is interacting with the dense envelope. At a higher resolution of 20 mas, or 44 au at the distance of G17, the continuum observations of Maud et al. (2019) uncovered a ring-like structure in a 120 au radius disk. They also found Keplerian kinematics seen in highly excited vibrational lines of water, tracing the hot upper layers of the disk.

G339.88-1.26. In their observations of the MYSO G339.88-1.26, Zhang et al. (2019a) found that within the midplane, CH₃OH and H₂CO trace the infalling, rotating envelope as well as the centrifugal barrier at 530 au, whereas SO₂ and H₂S trace the centrifugal barrier and outer disk. The radial difference between these two groups of molecules can be attributed to a radial temperature dependence and higher upper energy levels for the molecular transitions observed for the second group. In comparison, they found that SiO traces the Keplerian disk as well as the envelope and jet. They suggest that the enhancements in emission close to the centrifugal barrier are either explained by an accretion shock or the irradiated inner edge of the envelope. Outside of the midplane, CH₃OH emission is also found to be coincident with methanol masers. The thermal and maser methanol emission both extend perpendicular to the disk midplane in the outflow direction, lying to the south of the source on the blueshifted side of the disk (see their Figure 12(a)). Zhang et al. (2019a) explain this morphology, which is very similar to that seen for AFGL 4176, by shocks along the cavity walls of the outflow.

The blue-dominant lines in AFGL 4176, which include CH₃OH and H₂CO and have a ring-like morphology that is brighter on the blueshifted side of the disk, may be explained by similar processes to those seen in G339.88–1.26. Csengeri et al. (2018) also found that CH₃OH is tracing the centrifugal barrier in the MYSO G328.2551–0.5321. Thus, the enhancement of CH₃OH (and other blue-dominant molecules) at a specific radius from AFGL 4176 may be due to accretion shocks in the centrifugal barrier. In addition, the blueshifted asymmetry in the blue-dominant molecules may be either due to self-absorption of the redshifted emission from the centrifugal barrier by the envelope or by an actual asymmetry in the accreting material. Given that no blueshifted asymmetries are seen in our model for the lower K transitions of CH₃CN (J15), which would include the effect of self-absorption by the envelope, and that maser emission is seen coincident with the blueshifted CH₃OH emission (as noted also by Zhang et al. 2019a for G339.88–1.26), we suggest that there is an intrinsic asymmetry in the AFGL 4176 envelope and disk structure. This is also corroborated by the asymmetric NH₃(5,5) emission, which also peaks on the eastern side of the disk.

Looking back on the findings to-date for MYSOs with disks summarized above, it is clear that patterns are beginning to emerge, including the ring-like emission such as SO₂ and CH₃OH tracing shocks, possibly at the disk's centrifugal barrier and that N-bearing species often trace the disk. Nevertheless, there is still a large degree of inhomogeneity in the chemistry of these disks and their inner envelopes, and thus, observations of more sources at 100–1000 au scales are necessary to determine if any patterns persist and whether others emerge.