ALMA View of the Infalling Envelope around a Massive Protostar in S255IR SMA1

In Figure 1 we display in false-color the 900 µm dust continuum emission, which has a rms noise level of $\sim$ 1 mJy Beam^-1 and exhibits its strongest compact emission peaking at R.A. = 06:12:54.013 and Dec. = 17:59:23.05, the nominal position of the S255IR SMA1 source. The lower level continuum emission displays an elongation in the general northwest-southeast direction but also with a protrusion pointed towards northeast. Given the dust and gas are likely closely coupled and the gas temperature is at a level above 300 K toward the peak (see Section 3.3), the peak continuum intensity of 0.237 Jy Beam^-1, corresponding“ to a brightness temperature of 121 K, implies that the dust continuum emission, though with a non-negligible opacity, is mostly optically thin at this 0$\farcs$14 scale. With the assumption of a dust temperature of 300 K (see Section 3.3), we translate this peak intensity $F^{\mbox{beam}}_{\nu}$ to a gas column density $N_{\mbox{\tiny H${}_{2}$}}$ of 9.75 $\times$ 10²⁴ cm^-2 based on Equation A.26 in Kauffmann et al. (2008) as shown below

where $\Omega_{\mbox{\tiny A}}$ is the beam size of 0$\farcs$14, $\mu_{\mbox{\tiny{H2}}}$ is the molecular weight per hydrogen molecule of 2.8, $m_{\mbox{\tiny{H}}}$ is the atomic mass of hydrogen, $\kappa_{\nu}$ is the specific absorption coefficient (per mass), and B_ν(T) is the Planck function at a temperature T of 300K. Here, we adopt $\kappa_{\nu}$ of 0.010 cm²g^-1 by using $\kappa_{\nu}=\kappa_{250\mu m}({250\mu m}/{\lambda})^{\beta}$ with $\kappa_{250\mu m}=0.1$ cm²g^-1 (Hildebrand, 1983), $\beta=1.8$ (Chen et al., 2016), and implicitly a gas-to-dust mass ratio of 100. When the source distance of 1.78 kpc is taken into account, it implies enclosed within the central 250 au a total mass of around 0.36 M_⊙, or equivalently an average gas density of $\sim$ 6.6 $\times$ 10⁹ cm^-3. The integrated fluxes within larger circular apertures with diameters of 0$\farcs$5 ($\sim$ 900 au), 1$\farcs$0 ($\sim$ 1800 au), and 1$\farcs$5 ($\sim$ 2700 au) are 0.613 Jy, 0.932 Jy, and 1.08 Jy, respectively. With the same assumptions above, the enclose masses within thousand au scale are estimated to be 0.930 M_⊙, 1.41 M_⊙, and 1.64 M_⊙. These enclosed masses follow roughly $r^{0.5}$, where $r$ is the radius enclosing the volume, hinting a volume density profile scales roughly with $r^{-2.5}$ if a simple spherical geometry is considered. If a flattened two-dimensional (disc/toroid-like) geometry is assumed, the surface density could instead scale as $r^{-1.5}$ and the corresponding volume density in the inner region likely would be even higher than that quoted above due to the assumed geometry. We note, though, that the averaged temperature over a larger area would likely be cooler, and possibly certain missing flux may start kicking in, particularly when the angular scale reaches beyond 2″, the maximum recoverable angular scale of the observations. Both effects may lead to larger masses at larger volumes and the underlying density profile could hence be somewhat shallower than that inferred above. The overall mass estimates, nevertheless, suffer from various uncertainties as discussed in Section 4.4.

Figure 2 shows the spectrum spanning over the full spectral window within a synthesized beam centered at the continuum peak. We identified the most probable carriers of the observed spectral features in both our spectrum and image-cube by cross-referencing Splatalogue ¹¹1https://splatalogue.online/, an online database which consolidates the Jet Propulsion Laboratory Molecular Spectroscopy (JPL²²2https://spec.jpl.nasa.gov/, Pickett et al., 1998) and the Cologne Database for Molecular Spectroscopy (CDMS³³3https://cdms.astro.uni-koeln.de/, Müller et al., 2005) catalogs, and list them (in the order of descending frequency) in Table 1. The majority of these identified features originate from known interstellar complex organic molecules, including NH₂CHO, CH₃CHO, CH₃OCHO, CH₃CH₂OH, and (CH₂OH)₂, with multiple transitions at reasonable excitation energy levels present within the band, although analyzing their excitation is beyond the scope of this work. The remaining, generally weak, spectral features from our best guess but uncertain carriers are marked with question marks in the Table. The CH₃OH maser newly identified by Zinchenko et al. (2017) is clear seen as the strongest feature in the spectrum, with its intensity over 4 Jy Beam^-1(or equivalently over 2000 K in brightness temperature). This maser emission, at the strong peak location at the north-east of the continuum peak, as noted in Zinchenko et al. (2017), reaches $\sim$ 5.9 Jy Beam^-1at $\sim$ 0$\farcs$13 resolution and exceeds equivalently a brightness temperature of 3900 K (Zinchenko et al., 2017).

In this work, our primary focus is the CH₃CN transitions with their quantum numbers running from $J$=19–18, $K$=0 up to $K$=10, which are labeled in Figure 2. Among the detected molecular lines listed in Table 1, we highlight with boldface these transitions, whose upper level energies range from 167.7 K to 880.7 K. To better visualize the line profiles of the CH₃CN $J$=19–18 $K$ components in detail, we present in Figure 3 the individual $K$=0 to $K$=10 spectra extracted from Figure 2 but plotted in the velocity space. The typical observed full-width-half-maximum (FWHM) linewidths of the CH₃CN features are on the order of 8 km s^-1. As one can notice from both Figure 2 and Figure 3, the CH₃CN $K$=0 and $K$=1 transitions are blended themselves. A few other CH₃CN lines with higher $K$ are also noticeably blended by interlopers. Although the identifications of blended lines may not be obvious in this display (Figure 2), they become clearer in the position-frequency plot discussed later. The CH₃CN $K$=5 and $K$=6 transitions, for example, are blended with the CH${}_{3}^{13}$CN transitions, respectively, on the red-shifted and blue-shifted ends. Despite the line blending, one can spot in most CH₃CN transitions doubly peaked, asymmetric, (blue-)skewed profiles. The $K$=8 line is the only exception, where its red-shift end is contaminated by a relatively strong NH₂CHO feature. Doubly peaked blue-skewed profiles, if not resulted from line-blending, could be due to two velocity components, or as are often considered, an indication of inflow motion.

With a quick glance over the spectra, one may also find the fact that the intensities of the $K$=3 and even $K$=6 transitions are quite comparable to those of the $K$=0, 1, and 2 components. Given that the transitions of CH₃CN with $K=3n$ ($n=1,2,...$) are doubly degenerate as compared to the $K\neq 3n$ lines, the $K$=3 and $K$=6 components, for example, in the optically thin regime would presumably appear roughly twice brighter than their neighboring K lines. The comparable intensities in the observed spectra thus indicate that these lines are optically thick. A parallel line of evidence for the high optical depth of the lower $K$ lines comes from the observed CH${}_{3}^{13}$CN intensities as compared to those of CH₃CN. Although the isotopic ratio between ¹²C and ¹³C appears to vary with Galactocentric radii with large intrinsic scatters, the value is likely no smaller than 50 (Wilson & Rood, 1994; Wilson, 1999) for the S255IR region. Therefore, the column densities of CH₃CN and CH${}_{3}^{13}$CN and the corresponding peak intensities of their transitions, in the optically thin regime, should display a similar ratio. The intensity ratios of the observed CH₃CN $J$=19–18 $K$=0 and $K$=1 lines to those of the ¹³C isotopologues, as can be estimated from Fig 2, are much smaller than that, hinting again a significant opacity in the CH₃CN low $K$ lines. With the continuum baseline taken into account, the low-lying optically thick CH₃CN transitions appear saturated at an apparent brightness temperature of around 274 K at $\sim$ 0$\farcs$14 resolution.

The general distribution of the CH₃CN emission is depicted in the integrated intensity map of the CH₃CN $J$=19–18 $K$=3 transition shown by the contours in Figure 1. This integrated emission, as shown in the figure, is resolved with its strongest peak closely coinciding with the continuum peak. The map exhibits a prominent elongation of CH₃CN emission, which has an extent of around 1″ (or 1800 au) and a position angle (P.A., measured from north towards east) of 165 degree. This P.A. is nearly perpendicular to the axis of the large scale bipolar molecular outflows seen by the SMA and the IRAM 30M telescope with a P.A. = 67 degree (Zinchenko et al., 2015) as marked by the blue and red arrows in Figure 1, and is closely orthogonal to the small scale H₂O maser jet structure (with a P.A. = 49 degree) imaged through the Very Long Baseline Interferometric (VLBI) observations by Goddi et al. (2007) and Burns et al. (2016) using the Very Long Baseline Array (VLBA) and the VLBI Exploration of Radio Astrometry (VERA), respectively.

The full CH₃CN emission morphology, however, is more complicated. Similar to the continuum emission, there exists weaker level of CH₃CN emission extending toward the northeast direction. The channel maps of the CH₃CN $J$= 19–18, $K$=3 transition, shown in Figure 4, further detail the distribution of the CH₃CN emission around S255IR SMA1 at different velocities. From around 3 km s^-1 to 9 km s^-1, extended CH₃CN appears to be associated with the continuum and also toward the north-east. Part of the emission on the north-east end at around 3 km s^-1 and 5 km s^-1 may be co-spatial with the bright masing CH₃OH emission, which resides in this north-east region. At the extreme velocities ($<0$ km s^-1 or $>10$ km s^-1), the CH₃CN emission appears very compact in the close vicinity of SMA1.

Since CH₃CN is a symmetric-top rotor, radiative transitions between different $K$ levels are prohibited, and these levels are therefore populated solely by collisional processes. For this reason, the line ratios of these $K$ components are commonly used for probing the gas (kinetic) temperature. We perform synthetic spectral profile fitting of the CH₃CN spectra pixel by pixel following the approach in Wang et al. (2010) and Hung et al. (2019). One additional consideration for the analysis in this study is to take into account the contribution from the non-negligible continuum baseline in the true line brightness temperature. In the process we first carry out a simultaneous multiple–Gaussian profile fitting to the set of observed CH₃CN transitions based on the known frequency separation of these CH₃CN components, through which the systemic velocity and linewidth of the CH₃CN emission get determined. The process further incorporates a set of parameters, including the CH₃CN column density, gas (kinetic) temperature (which is assumed to be the same as the dust temperature), beam filling factor, and (dust) continuum optical depth for generating the synthetic spectra. The local thermal dynamical equilibrium (LTE) condition is assumed for setting the populations of different $K$ levels. By varying the input parameters and minimizing the difference ($\chi^{2}$) between the observed and the synthetic spectral profiles, we found the optimal set of parameters. The parameter space we considered is the following: 10¹⁴ cm^-2 – 10¹⁹ cm^-2 for the CH₃CN column density, 3 km s^-1 – 6 km s^-1 for the LSR velocity, 3 km s^-1 – 6 km s^-1 for the linewidth, 10 K – 500 K for the gas temperature, 0.02 – 1.0 for the filling factor, and 0.0 – 1.0 for the (dust) continuum optical depth. We limited our fitting exercise to the pixels where the CH₃CN $K$=3 transition is detected above 3$\sigma$ noise level.

Figure 5 presents the parameters inferred from the synthetic spectral profile fitting. Panel (a) shows the dust continuum brightness temperature, which is in essence equivalent to the continuum emission map shown in Fig. 1. Panel (d) displays the distribution of the CH₃CN column density. The dominant feature is a high column density ridge elongated in the northwest-southeast direction surrounding SMA1. From panel (e), we find that the (CH₃CN) gas temperature peaks toward SMA1 at more than 400 K and remains hot in the high column density regions. While the temperature gradually falls off toward the peripheral areas, within most of the fitted regions with an extent of 1$\farcs$ or 1800 au the temperature is greater than 100 K.

Based on the continuum emission brightness and the inferred gas/dust kinetic temperature, the dust optical depth , $\tau_{\mbox{\scriptsize dust}}$, is constrained and shown in panel (b). A clump with elevated opacity is apparent and coinciding with SMA1, with the dust continuum emission being mostly optically thin throughout the region when averaged over the synthesized beam. The inferred dust optical depth, based on equation A.9 and A.10 in Kauffmann et al. (2008), can be translated into the molecular hydrogen column density. We present the column density map in panel (c), which is fundamentally equivalent to panel (b) but scaled into the (log) column density unit. Toward the SMA1 peak at 250 au scale, we reach a molecular gas column density of $1.4\times 10^{25}$ cm^-2, consistent with the column density derived in Section 3.1 within 50%. As noted in Section 3.1, the derived mean gas density is on the order of 10⁹ – 10¹⁰ cm^-3. For the CH₃CN $J$=19–18 transitions, their critical densities for thermalization are around 10⁷ cm^-3. The “effective excitation densities” could be further reduced in the case of optically thick emission (Shirley, 2015). The LTE condition used in the spectral profile fitting is therefore a valid self-consistent assumption.

We attempt to establish empirically the radial profiles of the molecular gas column density and the (CH₃CN) gas temperature by extracting their values along the major axis from Figure 5 panels (b) and (e) and plotting them in Figure 6. Both profiles appear resolved as compared to the beam profile also shown in the same figure. We assume the gas column density, or equivalently the gas surface density, varies as $r^{-p}$ where $r$ is the radius. We find that the surface density profile, when convolved with the beam profile, is best fitted with $p=1.1$, as shown in the left panel of Figure 6. For the gas temperature, we adopt the CH₃CN temperature, which appears plateaued near the central YSO but vary with $r^{-0.7}$ beyond around 200 au. The (beam-convolved) fitted temperature profiles is plotted in Figure 6 right panel.

The inferred CH₃CN column density toward the SMA1 is $3.0\times 10^{17}$ cm^-2, thus implying a CH₃CN fractional abundance of $2\times 10^{-8}$. The same level of fractional abundance is found throughout the region as shown in Figure 5 panel (h). The value is somewhat higher when compared with the CH₃CN abundance, $10^{-9}$ (Su et al., 2009, and references therein), often found in massive star forming regions. As Su et al. (2009) revealed in the molecular envelope of the HII region G5.89, an enhancement of CH₃CN fractional abundance higher than $10^{-9}$ can occur in the inner and hotter region due to the evaporation of grain icy mantles. Indeed, this abundance is perfectly consistent with values found towards a sample of 17 hot molecular cores by Hernández-Hernández et al. (2014).

The CH₃CN velocity and linewidth (velocity dispersion) signatures are shown in Figure 5 panels (f) and (g), respectively. A velocity gradient is evident along a stripe from the northwest in red to the southeast in blue. To illustrate this velocity gradient across SMA1, we present in Figure 7 the position-frequency (P-F) diagram of spectral features within the entire spectral window made with a cut running through the SMA1 peak at a P.A. of 165 degree. The spectrum toward the continuum peak shown in Figure 2 is chiefly the intensity profile along the zero offset pixel in the cut. It is evident that, as those of CH₃CN (marked with white lines), all spectral features manifest very similar patterns, which enables better line identification by rejecting random noisy features along just one pixel. The carrier of these spectral features identified from this figure, including CH${}_{3}^{13}$CN and other species, are marked in yellow and red, respectively and also listed in Table 1. The few features with uncertain carriers noted earlier are further marked in orange. Still, there remains some features whose carriers are not identified.

A similar gradient in the general northwest-southeast direction has been reported by Wang et al. (2011) and Zinchenko et al. (2015). Wang et al. (2011) detected a velocity gradient in C¹⁸O and CH₃OCHO. The extend of the C¹⁸O emission, spanning nearly 10″ scale, was suggested to trace a rotating toroid that got fragmented into a multiple system including at least SMA1 and SMA2. The CH₃OCHO emission, not well resolved by their 1$\farcs$8 beam, exhibits a linear velocity gradient within 2″ scale and was thought to show likely the presence of a rotating and infalling core similar to the toroids described in Cesaroni et al. (2007). Zinchenko et al. (2015) interpreted the velocity gradient, observed at a higher resolution of 0$\farcs$5 in CH₃OH and CH₃CN, as a hint of Keplerian rotation within a rotationally-supported disc. In Figure 5 panel (f), we notice a hint of twisting in the position angle of this velocity gradient from the outside connecting toward the inner region around SMA1. The region with the highest velocity dispersion, as illustrated by Figure 5 panel (g), is very close to SMA1. In addition, a red-shifted component is visible toward the northeast, which is probably related to the red-shifted outflow lobe.

From Figure 7, we extract and plot in Figure 8 the position-velocity (P-V) diagram of the CH₃CN $J$=19-18, $K$=0 through $K$=10 transitions for further scrutinizing the northwest-southeast velocity gradient across SMA1. In the first panel, the $K$=0 and $K$=1 components are blended and shown together. Apparent blending features are also discernible in the $K$=4, $K$=8, and $K$=9 panels. Two important attributes in the P-V diagram are readily recognized. First, an overall velocity gradient across the cut is again evident from the blue-shifted gas at the positive (SE) end to the red-shifted gas at the negative (NW) end. There is a quite abrupt flip of the velocity near the central position. Second, while the blended $K$=0 and $K$=1 components appear relatively extended, tracing a smaller velocity gradient, at progressively higher $K$s, the emission extent becomes more compact and the apparent velocity gradient gets steeper.

What is the origin of this observed velocity gradient pattern? Since this gradient is closely perpendicular to the general direction of the bipolar outflows, it most plausibly indicates rotation. Furthermore, as the higher K transitions trace the more compact, hotter, and presumably denser gas closer to the embedded protostellar object, the increasing velocity gradient signifies a spin-up rotation toward the center. It is tempting to attribute the velocity gradient to the Keplerian rotation of a rotationally supported circumstellar disc around the central YSO. Plotted in Figure 9 panel (a) are, respectively, red, blue, and cyan curves corresponding to Keplerian rotation of (8.04/$sin^{2}i$) M_⊙, (3.57/$sin^{2}i$) M_⊙, and (0.89/$sin^{2}i$) M_⊙, where $i$ represents the inclination angle of the system with $i=0$ meaning face-on. For an inclination angle of 25 deg, the red, blue, and cyan line, for example, denote the Keplerian curves for a central mass of 45 M_⊙, 20 M_⊙, and 5 M_⊙, respectively. For an inclination angle of 60 deg, the red, blue, and cyan line, then stand for the Keplerian curves for a central mass of $\sim$ 10.7 M_⊙, 4.8 M_⊙, and 1.2 M_⊙, respectively. We plot in Fig 1 the perimeters of a axisymmetric circumstellar structure with a diameter of 1″, a P.A. of 165 degree, and an inclination angle of 0 degree, 25 degree, and 60 degree. The degree of elongation in our integrated CH₃CN emission, namely the minor axis of the structure being roughly half in length of the major axis, would suggest an inclination angle of no less than 60 degree. It is evident that the emission pattern on the red-shifted and blue-shifted sides are not fully symmetric. Meanwhile, the blue curve matches reasonably the emission pattern at the two ends. When viewed at an inclination of 60 degree, this curve, as noted above, would correspond to the Keplerian motion around a central mass of $\sim$ 4.8M_⊙.

To aid the interpretation of observed P-V signatures, we pursue a set of radiative transfer model simulations using the ”Simulation Package for Astrophysical Radiative Trans(X)fer” (SPARX) ⁴⁴4https://sparx.tiara.sinica.edu.tw/. The software package produces image cubes through radiative transfer calculations of molecular line radiation with input physical models. We further convolve these synthetic image cubes output by SPARX with the observing resolution (0$\farcs$14) using the ”imsmooth” command in CASA and generate the model P-V diagrams with the CASA ”impv” task. We note that no interferometric spatial filtering was applied to form the model cubes during this process as we do not expect significant missing flux within the angular scale ($\sim$ 1″) of the observed structure.

For the physical model setup, we first assume that the emission originates from a flattened two-dimensional structure, possibly a disc or an envelope. This is consistent with the elongation we see in the continuum as well as in the CH₃CN integrated emission. For the molecular gas column density and temperature, we adopt the empirical profiles inferred in Section 3.3. The column density radial profile is scaled to match the the observed gas column density or equivalently the observed (dust) continuum intensity with the gas temperature and a gas-to-dust ratio of 100 factored in. We then set a CH₃CN fractional abundance be $2\times 10^{-8}$ as suggested by Figure 5 and discussed earlier. The full size of the system is not fully known but we truncate the CH₃CN abundance at a radius of 900 au ($\sim 0\farcs 5$) based on the extent of the CH₃CN emission observed in the P-V diagram. We carry out our calculations in the LTE regime since the observed transitions are most likely thermalized.

In Figure 9 we compare the observed $J=19-18,K=3$ emission features (shown in false color in all panels and in black contours in panel (b)) with the modeling results. In toy Model 1 (M1) and Model 2 (M2), we assumed the emissions are from gas following pure Keplerian rotation. The blue contours in panel (c) for M1 and panel (d) for M2 display the P-V patterns with their central YSO masses and inclination angles being 20 M_⊙and 4.76 M_⊙, and 25 degree and 60 degree, respectively. As one can see, the outer most contours of the model emission traces reasonably well the extent of the observed emission at the same level, suggesting that the empirical density, temperature, and abundance values helped gauging the simulation intensity be in the right ballpark. There are slight differences between the two model patterns at small radii, but they generally agree with each other and resemble the observed features. Nevertheless, there are subtle sub-structures that are not well reproduced. In particular, the central emission ”waist” between the two ends at the two ”quadrants” in the simulated P-V pattern is significant narrower than that of the observed emission. The observed wider ”waist” features can be understood as signatures of infalling motions, which have been modeled and demonstrated extensively for low mass YSO systems by, for example, Ohashi et al. (1997); Yen et al. (2013); Sakai et al. (2014a). Such infalling features have also been recently outlined at high angular resolution toward, for example, the massive YSO G023.01–00.41 (see Fig. 5 of Sanna et al. (2019)). In panel (e), we subsequently explore an infalling-rotation envelope model in which the gas moves in the azimuthal direction following the Keplerian manner (around a central mass of 4.76 M_⊙at an inclination angle of 60 degree) but also bears a radial infall velocity of 2.5 km s^-1at all radii. The P-V emission pattern from this envelope model (M3), shown in Figure 9 panel (e), captures the observed features well. In particular, the infall motion leads to a wider ”waist” in the P-V diagram. Meanwhile, as opposed the M1 and M2 cases in which the emission peaks in the ”central spine” of the ”waist”, the emission is stronger on the ”edges” of the ”waist”, similar to that in the observed pattern. In Figure 9 panel (f), we furthermore hypothesize another model (M4) in which the gas follows the accretion flow prescribed by Ulrich (1976). The essence of such an accretion flow, which has been adopted for studying the envelopes of low- and high-mass stars forming regions in, for example, Yen et al. (2014) and Keto, & Zhang (2010), is that the gas parcel, while accelerates toward the central YSO due to the gravity, also have a rotational motion along a (ballistic) trajectory with its angular momentum conserved. In this scenario, the gas stream lines are dictated by two of the three parameters, the (specific) angular momentum, the central mass, and the so called “centrifugal barrier radius”, at which point the infalling material can not move further in but presumably would get shocked, lose angular momentum, and settle into Keplerian rotation orbits around the central mass. Admittedly, we have no knowledge for this model setup about the (specific) angular momentum or the centrifugal barrier radius, which is, however, likely to be small and unresolved. For the purposing of testing, we assume arbitrarily its centrifugal barrier at a small radius of 90 au (or 0$\farcs$05). As the gas density and thermal structures are readily parameterized with the empirical profiles, we mainly adopt the velocity fields of the accretion flow. In this case, as in model M3, the feathers matches those with the observations quite well - the emission pattern is wider at the ”waist” and is stronger on the edge of that. The parameters of all the four radiative transfer models are shown in Table 2.

In brief, models M3 and M4, with infalling gas motion factored in, exhibit better correspondence to the features seen in the P-V diagram as compared to models M1 and M2, in which only Keplerian motion is considered. Given the reasonable match of M4, we present in Figure 10 the P-V diagram of the observed CH₃CN emission for all $K$ components in false color and the corresponding model emission in contours. These contours can be directly compared with those in Figure 8 and their consistencies can be recognized.

To interpret the observed CH₃CN kinematic signature, it is imperative to have some knowledge about the mass of the central YSO NIRS3 and the geometry of the system. Based on a distance of 2.5 kpc and the SED analysis, Ojha et al. (2011) previously found a mass of $\sim$ 27 M_⊙ for the central exciting star in SMA1. By adopting an updated distance (of $\sim$ 1.6 kpc) hence a lower luminosity, Zinchenko et al. (2015) revised the mass of the central YSO in SMA1 to be $\sim$ 20 M_⊙. Indeed, considering the theoretical calculations by, for example, Hosokawa et al. (2016) in which both stellar luminosity and accretion luminosity (with a mass accretion rate of $10^{-3}$ to 10^-4 M_⊙) are included, one would translate the observed luminosity of S255IR NIRS3/SMA1 at a level of at least $2.4\times 10^{4}$ L_⊙ (e.g. Caratti o Garatti et al., 2017) to a mass of $\sim$ 15 – 20 M_⊙ for the central YSO.

While in Section 3 we have adopted an inclination angle of 60 degree for the radiative transfer models M2, M3 and M4, this angle in fact has not been very well constrained observationally. According to the velocity gradient seen in the CH₃OH emission, Zinchenko et al. (2015) previously inferred an inclination angle of $\sim$ 25 degree for an Keplerian disc surrounding the central 20 M_⊙YSO in SMA1. The elongated morphology of the CH₃CN emission seen by our observation, however, does not fully support such a scenario but instead would imply an inclination angle no less than 60 degree. Incidentally, the mid-infrared (N band) observations with the MIDI instrument on the Very Large Telescope Interferometer (VLTI) carried out by Boley et al. (2013) implied a highly inclined (79.6 degree) disc-like structure with an extension of 92.8 mas (or equivalently 165 au) in FWHM oriented at a P.A. of 149 degree based on their visibility fitting calculation with a two-dimensional Gaussian model. The finite angular sampling and coverage in the $u$-$v$ plane (only between 32 degree and 76 degree) of their observations, however, might have limited the reliability of the measurement. Besides, Burns et al. (2016), based on the maser kinematics measured both along the line-of-sight and in the sky plane through VLBI observations, suggested a jet from SMA1 inclined at 86 degree with a P.A. of 49 degree. This jet geometry would imply also a highly inclined (86 degree) disc with a P.A. of around 140 degree, similar to that inferred from the VLTI MIDI observations. In short, the orientation of the major axes obtained from the majority observations, including ours, are generally consistent. Meanwhile, although the circumstellar disc has not been directly imaged, the inferred angle of such a structure is likely highly inclined. The value of 60 degree assumed by us appears to be a lower limit. Since the YSO mass is reasonably between 15 – 20 M_⊙, and the rotation axis is inclined by more than 60 degrees with respect to the line-of-sight (as also suggested by the well collimated outflow driven by the central star), the assumption of Keplerian rotation does not reproduce the observations properly and has to be revised.

Numerical radiation hydrodynamical simulations of massive star formation such as the one introduced in the above section or those carried out in, for example, Klassen et al. (2016) and Meyer et al. (2017) suggest that the formation of disc, as it evolves, may appear highly asymmetric with clumps and/or spiral forms due to gravitational instability at hundreds to even thousands of au scales. Gas clumps at these scales subsequently lose angular momentum and fall in towards 10s au scales and finally accrete on to the central massive YSO, giving rise to markedly elevated luminosity. While we do not witness direct signatures of spiral arms, our CH₃CN data do suggest inhomogeneity and asymmetric kinematic features in the envelope, such as the asymmetric emission pattern seen in the P-V diagram. The limited physical resolution of $\sim$ 250 au prevented us from revolving and identifying the small (10s au) scale gas clumps that might have been physically associated with or left over from the most recent burst event.

Johnstone et al. (2013) computed the temperature and luminosity variation of deeply embedded (low mass) stars due to episodic accretion events. They demonstrated that the dust temperature responses much quicker than the gas temperature to luminosity variations. In particular, at the outer envelope where the dust emission is optically thin, heating of the gas would take year-long timescale as compared to that (month-long timescale) for the dust. For S255IR SMA1, evidences show that the burst occurred in 2015 has been waning down (Liu et al., 2018; Szymczak et al., 2018; Uchiyama et al., 2019). While our CH₃CN data from the 2017 Sep (post burst) epoch do not show changes in their kinematic signatures, there appear tentative hints that the line intensities may have tailed off, suggesting reduced gas temperature in the region. High resolution follow-up observations of CH₃CN in the post-burst phase shall better examine this possible signature of gas temperature variation.

We have estimated the molecular column density and gas mass toward the S255IR SMA1 by using the dust continuum emission. There exist a few potential caveats in inferring the gas mass due to the uncertainties in the assumed dust temperature, opacity, and measured flux. As shown in Liu et al. (2018), the continuum flux of this region was varying during the luminosity burst. The changing continuum flux is presumably reflecting the changing temperature, or even possibly dust properties including its opacity. Only when the flux, temperature, and dust opacity are secured will we able to confidently measure the column density. In the case of the massive star forming region G351.77-0.54, the gas temperature derived from the CH₃CN $J=37-36$ $K$-ladder was different from the dust temperature inferred from the dust continuum (Beuther et al., 2017). It was suggested by the authors that the dust continuum emission may trace the deeper and colder region while the CH₃CN emission originates from the hotter surface layer. In our case, it is also possible that the gas and dust temperatures get de-coupled during a burst event as suggested by Johnstone et al. (2013). In addition, we have adopted the dust absorption coefficient $\kappa_{\nu}$ of 0.010 cm²g^-1 based on Hildebrand (1983) and $\beta=1.8$ from Chen et al. (2016). The index $\beta$ is, however, highly uncertain — a change of $\beta$ to a smaller value of 1.5 or a larger value of 2 would, respectively, lead to an increase of $\kappa_{\nu}$ by 46% or a decrease of $\kappa_{\nu}$ by 23%. Alternatively, if we consider the value of dust absorption coefficient for ”naked” dust at a density of 10⁸ cm^-3 as tabulated in Table 1 of Ossenkopf & Henning (1994) and again a gas-to-dust ratio of 100, $\kappa_{\nu}$ would increase by over a factor of 10 to be 0.113 cm²g^-1, thus reducing the mass of the infalling envelope and increasing the CH₃CN abundance within by the same factor. Finally, the observed flux variation in the 850 $\mu m$ continuum in the luminosity burst is about a factor of 2 (Liu et al., 2018). Considering all these aspects, the uncertainly in the mass estimate could be a factor of a few. The currently derived gas mass of around 1.64 M_⊙ within $\sim$ 2700 au appears significantly less than the presumed central stellar mass of $\sim$ 20 M_⊙, and hence the subsequent model calculations of gas kinematics without considering their self-gravity should be reasonable.

Chemical differentiation in between the circumstellar disc and envelope components has been seen toward low-mass YSOs, and it is perhaps best exemplified by the case of the Class 0 YSO L1527. Sakai et al. (2014a, b) found that toward L1527 the (cyclic-)C₃H₂, CCH and CS emissions are doubly peaked and tracing the infalling rotating envelope. On the other hand, the SO, H₂CO, and CH₃OH emissions are centrally peaked, tracing primarily the rotationally supported Keplerian disc component inside the centrifugal barrier. Lee et al. (2019), by vertically resolving the almost edge-on protostellar disc in the low-mass YSO HH212, reported that complex organic molecules such as CH₃OH, CH₃CHO, CH₃CH₂OH, as well as H₂CCO, HCOOH, and CH₃OCHO, normally considered as tracers of hot molecular cores and corinos, seem to reside in the warm disc atmosphere within the centrifugal barrier.

Similar signatures of chemical differentiation have also been indicated toward massive YSOs. For example, in the study of G339.88-1.26 by Zhang et al. (2019), there appear distinct spatial extent and kinematic signatures in the position-velocity diagrams for different species. In this study, CH₃OH and H₂CO are tracing a more extended infalling rotating envelope while SO₂ and H₂S are most likely tracing a disc inside and particularly enhanced at the centrifugal barrier (Zhang et al., 2019).

For the case of S255IR SMA1, we have detected and listed a good number of species/transitions in Table 1, as noted in Section 3.2. The majority features seem to all trace the same full extent of the envelope we see. Partially blended CCS and SO₂ emission features are detected at around 349.225 GHz within the band. While we do not find firm evidences of differentiation in terms of their patterns in the P-V diagram (see Fig. 7), there are hints that the CCS and SO₂ emissions may originate more prominently from smaller radii, which appears somewhat similar to the case of G339.88-1.26. Details of the full chemical signatures in this region will be further discussed in a forthcoming study, and observations with higher sensitivity are required to verify the case of chemical differentiation.