Discovery of a sub-Keplerian disk with jet around a 20Msun young star. ALMA observations of G023.01-00.41

It is well established that Solar-mass stars gain mass via disk accretion, until the mass reservoir of the disk is exhausted and dispersed, or condenses into planetesimals. Accretion disks are intimately coupled with mass ejection via polar cavities, in the form of jets and less collimated winds, which allow mass accretion through the disk by removing a substantial fraction of its angular momentum. Whether disk accretion is the mechanism leading to the formation of stars with much higher masses is still unclear. Here, we are able to build a comprehensive picture for the formation of an O-type star, by directly imaging a molecular disk which rotates and undergoes infall around the central star, and drives a molecular jet which arises from the inner disk regions. The accretion disk is truncated between 2000-3000au, it has a mass of about a tenth of the central star mass, and is infalling towards the central star at a high rate (6x10^-4 Msun/yr), as to build up a very massive object. These findings, obtained with the Atacama Large Millimeter/submillimeter Array at 700au resolution, provide observational proof that young massive stars can form via disk accretion much like Solar-mass stars.


Introduction
Models of massive star formation in the disk accretion scenario predict that circumstellar disks could reach radii between 1000 and 2000 au (e.g., Krumholz et al. 2007;Kuiper et al. 2011;Harries et al. 2017). Observationally, evidence for gas rotation near young stars with tens of Solar masses exists (e.g., Qiu et al. 2012;Johnston et al. 2015;Cesaroni et al. 2017), although the amount of gas mass undergoing rotation is a significant fraction of the star mass, and these envelopes, possibly hosting an inner disk, might be prone to fragmentation due to self-gravity and/or develop spiral instabilities (e.g., Kratter et al. 2010;Kuiper et al. 2011;Klassen et al. 2016;Chen et al. 2016;Meyer et al. 2017Meyer et al. , 2018. The interplay between disks and jets provides a mechanism to ensure mass accretion through the disk. Notwithstanding their connection, there is poor evidence of disk-jet systems in the inner few 1000 au of stars with tens of Solar masses (e.g., Beltrán & de Wit 2016), such as those resolved around Solar-and intermediate-mass stars (e.g., Lee et al. 2017a,b;Cesaroni et al. 2005Cesaroni et al. , 2013Cesaroni et al. , 2014. In order to establish the disk accretion scenario as a viable route for the formation of the most massive stars (e.g., Kuiper et al. 2010), we seek to simultaneously resolve the spatial morphology of the disk and the regions where the jet originates, and to determine whether or not the disk gas is in centrifugal equilibrium.
In the recent years, we have identified an "isolated" massive young star in the hot molecular core (HMC) of the starforming region G023.01−00.41, which is located about half way between the Sun and the Galactic center, at a parallax distance of 4.6 kpc (Brunthaler et al. 2009). The HMC emits a luminosity of 4 × 10 4 L (Sanna et al. 2014), and stands out among the strongest Galactic CH 3 OH maser sources (Menten 1991;Cyganowski et al. 2009;Sanna et al. 2010;Moscadelli et al. 2011;Sanna et al. 2015). The HMC luminosity corresponds to that of a zero-age-main-sequence (ZAMS) star with a mass of 20 M and an O9 spectral-type (Ekström et al. 2012). The HMC is located at the center of a collimated CO outflow whose emission extends up to the parsec scales (Araya et al. 2008;Furuya et al. 2008;Sanna et al. 2014); we tracked the driving source of the outflow down to the HMC center, where we imaged a collimated radio thermal jet associated with strong H 2 O maser shocks (Sanna et al. 2010(Sanna et al. , 2016. The different outflow tracers are aligned on the sky plane and constrain the direction of the outflow axis within an uncertainty of a few degrees (Sanna et al. 2016). An accurate knowledge of the outflow axis allows us to pinpoint the young stellar object (YSO) position, and thus to circumvent the usual problem of distinguishing between velocity gradients due to expanding or rotational motions. Notes. Column 1: granted 12m-array configuration. Columns 2 and 3: target phase center (ICRS system). Columns 4: source radial velocity. Columns 5: minimum and maximum rest frequencies covered. The ALMA IF system was tuned at a LO frequency of 226.544 GHz, and made use of 4 basebands evenly placed in the lower and upper sidebands. Column 6: maximum spectral resolution required. Columns 7, 8, and 9: bandpass, phase, and absolute flux calibrators employed at both runs. Calibration sources were set by the ALMA operators at the time of the observations. Columns 10: required beam size at a representative frequency of 220.6 GHz.
Here, we exploit the information about the star-outflow geometry in G023.01−00.41, and make use of Atacama Large Millimeter/submillimeter Array, ALMA, observations (Sect. 2), to directly resolve a disk-jet system associated with an O-type YSO. We firstly show that the line emission from dense gas reveals a molecular disk which extends up to radii of 2000-3000 au from the central star; the disk is warped in the outer regions (Sect. 3). Then, we make use of the position-velocity (pv-) diagrams of gas along the disk plane to show that gas is falling in close to free-fall, and slowly rotating with sub-Keplerian velocities; at radii near 500 au, gas rotation takes over, and approaches centrifugal equilibrium at shorter radii (Sect. 4). In the process, we image the molecular jet component which arises from the inner disk regions. We corroborate our conclusions by comparing the observed pv-diagrams with those obtained for a disk model around a 20 M (Appendix A).

Observations and calibration
We made use of our previous Submillimeter Array (SMA) observations at 1 mm (Sanna et al. 2014), which covered a range of angular resolutions between 3 and 0 . 7, to set up the requirements for higher resolution observations with ALMA.
We observed the star-forming region G023.01−00.41 with the 12 m-array of ALMA in band 6 (211-275 GHz). Observations were conducted under program 2015.1.00615.S during two runs, on 2016 September 5 and 16 (Cycle 3), with precipitable water vapor of 1.5 mm and 0.6 mm, respectively. The 12 m-array observed with 45 antennas covering a baseline range between 16 m and 3143 m, with the aim to achieve an angular resolution of 0 . 2, and to recover extended emission over a maximum scale of 2 .
We made use of the dual-sideband receiver, in dual polarization mode, to record 12 narrow spectral windows, each 234 MHz wide, and an additional wide band of 1875 MHz. The narrow spectral windows were correlated with 960 channels and Hanning smoothed by a factor of 2, achieving a velocity resolution of 0.7 km s −1 . Individual spectral windows were placed to cover a number of molecular lines of high-density tracers (> 10 6 cm −3 ) such as methanol (CH 3 OH) and methyl cyanide (CH 3 CN). These settings provide a line spectral sampling of 13 channels for an expected linewidth of 9 km s −1 . The wide band was correlated with 3840 channels and centered at a rest frequency of 217.860 GHz. This band was used to construct a pure continuum image from the line-free channels.
We spent 1 hour of time on-source during a total observing time of 2 hours, which includes calibration overheads. Time constraints were set to achieve a thermal rms per (narrow) spectral channel of about 2 mJy beam −1 (with 36 antennas), which corre-sponds to a brightness temperature of 1.4 K over a beam of 0 . 2. Additional observation information is summarized in Table 1.
The visibility data were calibrated with the Common Astronomy Software Applications (CASA) package, version 4.7.1 (r39339), making use of the pipeline calibration scripts. We determined the continuum level in each spectral window separately. We made use of the pipeline spectral image cubes and selected the line-free channels from a spectrum integrated over a circular box of 2 in size, which was centered on the target source. The task uvcontsub of CASA was used to subtract a constant continuum level across the spectral window (fitorder = 0). We imaged the line and continuum emission with the task clean of CASA. In each individual map, the ratio between the image rms and the thermal noise, expected from the ALMA sensitivity calculator, is near unity. For the continuum map, we integrated over a line-free bandwidth of 218 MHz selected from the wide spectral window, and achieved a signal-to-noise ratio of about 100. Imaging information is summarized in Table 2.

Results
According to the "observer's definition" of disk framed in Cesaroni et al. (2007), necessary conditions for a disk detection are (i) a flattened core of gas and dust perpendicular to the outflow direction, and (ii) a velocity gradient along the major axis of the core. We therefore want to analyze the spatial morphology and kinematics of dense gas which lies at the center of the bipolar outflow, and, specifically, along a position angle of −33 • on the sky plane (measured east of north). This direction traces the projection of the equatorial plane, hereafter the "disk plane", perpendicular to the molecular outflow and the radio jet orientation (or, simply, the outflow), which is inclined by less than 30 • with respect to the sky plane (Sanna et al. 2014(Sanna et al. , 2016. In order to gain an overall view of the mass distribution around the star, we first looked for a molecular tracer of cold and dense gas (> 10 6 cm −3 ). In the left panel of Fig. 1, we plot the velocity-integrated map of the CH 3 OH gas emission at a frequency of 218.440 GHz (in color); this CH 3 OH transition has an upper excitation-energy of only 45 K (E up ). To emphasize the bulk of the emission, we integrated over a small velocity range (2 km s −1 ) near the line peak.
For clarity, we have rotated this map in order to align the outflow axis, which lies at a position angle of +57 • on the sky plane (Sanna et al. 2016), with the vertical axis of the plot. The blueshifted (approaching) lobe of the outflow points now to the north. The horizontal axis is aligned with the disk plane, which is marked in the plot at multiples of 1000 au from the central star position. The position of the central star was set as to maximize the symmetry of the blue-and redshifted sides of the disk (see Figs. 2 and 4), and has coor- Left: momentzero map of the CH 3 OH (4 2,2 -3 1,2 ) E line emission (colors) combined with the continuum map of the dust emission at 1.37 mm (dashed white contours). Maps have been rotated clockwise by a position angle of −57 • , in order to align the (projected) outflow axis, drawn on the left side, with the north-south axis; negative outflow heights indicate the receding outflow direction. The CH 3 OH emission was integrated in the range 78.8-80.2 km s −1 ; the wedge on the top left corner quantifies the line intensity, from its peak to the maximum negative in the map. The lowest dashed contour corresponds to the 10 σ level of the dust map, and the inner contour traces the 50% level from the continuum peak emission. The disk plane is drawn at three radii: from the central star position (star) to 1000 au (yellow), from 1000 to 2000 au (black), and up to 3000 au (pink). The dotted black contour draws the 60% level of the CH 3 OH emission, which identifies the disk profile (see text). The synthesized ALMA beams, for the dust continuum map (dashed circle) and the line map (red circle), are shown in the bottom left corner. Right: first-moment map (colors) of the CH 3 OH (4 2,2 -3 1,2 ) E line emission plotted in the left panel. The outer contour traces the 40% level of the moment-zero map; the inner dotted contour is the same as in the left panel. The LSR velocity scale is drawn in the upper left.  Notes. Columns 1 and 2 list the tracer and central frequency of each map (or molecular transition), respectively. For the line emission, column 3 specifies the upper excitation-energy of the molecular transition. Column 4 reports the Briggs' robustness parameter set for the imaging. Column 5 reports the restoring (circular) beam size, set equal to the geometrical average of the major and minor axes of the dirty beam size. Columns 6 and 7 report the rms of the map and the peak brightness of the emission, respectively. a Units of mJy beam −1 km s −1 .
dinates of R.A. (J2000) = 18 h 34 m 40. s 283 and Dec. (J2000) = -9 • 00 38 .310 (± 30 mas per coordinate). This position also coincides with the peak position of the radio continuum emission observed with the Very Large Array at the same resolution as the ALMA beam (Sanna et al. 2016). The same geometry and symbols will be used in all maps for comparison. Figure 1 shows that dense gas condenses in the direction perpendicular to the outflow axis. The CH 3 OH gas emission extends up to a deconvolved radius of 2380 ± 840 au (or 520 mas) from the central star, and has a ratio between the minor-to-major axes of 0.2. These values are calculated by Gaussian fitting the CH 3 OH emission within the 60% contour level (dotted black contour in Fig. 1). We use this threshold to select the region where the CH 3 OH iso-contours have a nearly constant height. This piece of evidence satisfies the first condition, that of a highly flattened structure perpendicular to the outflow, and hereafter we refer to the 60% contour level as the "disk profile". Below the 60% contour level, the CH 3 OH emission emerges at (relatively) large radii from the star (>2000 au) and appears significantly warped with respect to the disk plane. In this paper, we focus on the analysis of the emission within the disk profile.
In the right panel of Fig. 1, we plot the intensity-averaged map of the velocity field of the CH 3 OH gas (i.e., the firstmoment map). There is a clear gradient between the redshifted and blueshifted velocities moving from the left to the right side of the star along the disk plane. This second piece of evidence fulfills the more stringent requirement for confirming a disk candidate. On a closer look, the velocity field is actually skewed along the main outflow direction. This is not surprising, since, around Solar-mass stars, CH 3 OH emission has been shown to be excited both in the disk and jet regions (e.g., Leurini et al. 2016;Lee et al. 2017c;Bianchi et al. 2017). Indeed, contamination of the velocity gradient across the disk plane, by the blueshifted and redshifted gas velocities due to the outflow expansion, to the north and south of the star respectively, can explain this apparent velocity gradient.
To better investigate the gas kinematics inside the disk profile, in Fig. 2 we plot a number of maps of the CH 3 OH (10 2,8 -9 3,7 ) A + line emission at different velocities. This transition has an excitation-energy (E up = 165 K) about four times higher than that of the line in Fig. 1. The upper row shows the spatial morphology of the CH 3 OH gas emission at increasing (redshifted) LSR velocities, from left to right, and underlines the eastern disk side. On the contrary, the lower row shows the CH 3 OH gas distribution at decreasing (blueshifted) LSR velocities, which arise from the western side of the disk. Each row spans a velocity range of 2.8 km s −1 , which is the rotation speed expected around a 20 M star at an outer radius of 2300 au (assuming centrifugal equilibrium). Figure 2 resolves the spatial distribution of gas around the star at the different velocities, providing direct proof that the CH 3 OH gas is rotating clockwise around the outflow axis (as seen from the north), with the approaching and receding sides of the disk to the west and east of the central star, respectively. The pair of maps on the same column to the left shows the transition between the eastern (upper) and western (lower) sides of the disk emission, and can be used to infer the rest velocity of the star (see below). Between velocities of 80 and 79 km s −1 , the CH 3 OH gas emission is maximally stretched along the disk plane on either sides of the star; at higher (lower) velocities, the emission progressively brightens southeastward (northwestward) of the disk plane, in agreement with the redshifted (blueshifted) outflow lobe. A posteriori, this behavior confirms that the observed molecular species can trace both the disk kinematics and the (inner) outflowing gas (see also Fig. 4, right).
In the Fig. 3, we provide the same analysis for the CH 3 OH (18 3,16 -17 4,13 ) A + line emission, which has an excitation-energy of 446 K (E up ), and is more affected by the outflow emission. For this line, the transition between the red-and blueshifted sides of the disk occurs at lower velocity than in Fig. 2. The rest velocity of the star, V , was estimated from the median value between Fig. 2 and Fig. 3, and is set to +79.1 km s −1 with an uncertainty of ± 0.4 km s −1 . We explicitly note that, since the disk plane is not seen exactly edge-on, the blueshifted outflow emission close to the star, which lies between the observer and the disk, shifts the velocity field towards the lower velocities, when weighting the velocity channels by the line intensity. In the first-moment map of Fig. 1, the net result is that the redshifted velocities are down-weighted, and the central velocity is shifted by about −1.5 km s −1 with respect to the rest velocity of the star.
Overall, in Figs. 1 and 2 we resolve the gas emission within a few 1000 au of the O-type YSO, which allows us to prove the presence of a disk, and prompts us to study its kinematics in detail.

Discussion
In the following, we want to study the dependence of the gas velocity on the distance to the star, and to quantify the disk mass. Together with the observations, we discuss the results of a radiation transfer model for a circumstellar disk around a 20 M star (Appendix A). We have simulated the disk appearance for our specific observational conditions, and produced synthetic maps of the line and dust emission around the young star, in order to compare the observed pv-diagrams with those expected under two simple assumptions: the disk is falling in towards the central star due to gravitational attraction; the disk is rotating around the central star in centrifugal equilibrium.

Position-velocity diagrams
In Fig. 4, we study the velocity profile of gas along the disk plane through the pv-diagrams of two molecular species having a common methyl group, CH 3 OH and CH 3 CN. On the same column to the left, we compare the pv-diagrams of the CH 3 OH (10 2,8 -9 3,7 ) A + and CH 3 CN (12 4 -11 4 ) line transitions, which have similar excitation-energies of 165 K and 183 K, respectively. The two plots span the same ranges in space and velocity, and the lower intensity contour is drawn as to include an outer radius of 2400 au. Both pv-diagrams are double peaked, and the peaks are symmetrically displaced with respect to the star position and velocity. However, the CH 3 CN line emission shows a steeper velocity gradient at small offsets from the star, which indicates that this transition traces a region closer to the disk center than the CH 3 OH transition. This piece of evidence also explains the lower redshifted tail of the CH 3 CN pv-diagram, which traces the higher rotation speeds approaching the central star, as it can be expected for Keplerian-like rotation (v rot ∝ R −1/2 ). This component is missing in the roundish contours of the CH 3 OH emission. In Fig. 5, we plot the pv-diagram of the CH 3 CN (12 3 -11 3 ) line, which has an upper excitation-energy of 133 K, and shows a similar profile to that of the K = 4 line. Their agreement, when compared to the pv-diagram of the CH 3 OH line, strengthens the idea that the CH 3 CN gas allows us to peer into the innermost disk regions, and this holds for a broad range of line excitationenergies. On the other hand, both the CH 3 CN and CH 3 OH lines show increasing blueshifted velocities moving close to the star. In the right column of Fig. 4, we show a channel map at the blueshifted velocity of +73.2 km s −1 for both transitions. These maps clearly show a compact (molecular) outflow component, namely the jet, to the north of the star, corresponding to the direction of the blueshifted outflow lobe. Since the blueshifted outflow emission lies in between the observer and the disk, we cannot neglect its influence on the pv-diagrams (see below). Theory predicts that jets are launched and collimated in the inner few 100 au from the central star (e.g., Frank et al. 2014;Kuiper et al. 2015Kuiper et al. , 2016. The jet width of 1510 ± 73 au (or 329 mas), determined from the deconvoled (Gaussian) size of the CH 3 OH emission, allows us to set an upper limit of 800 au to the radius where the jet originates. Notably, we do not detect the redshifted lobe of the jet, which lies in the background with respect to the disk, and interpret this result as evidence for the inner disk regions being (partially) optically thick at 1 mm (e.g., Sanna et al. 2014;Forgan et al. 2016).
In Fig. 6, we compare the observed pv-diagram of the CH 3 CN (12 4 -11 4 ) line, which is sensitive to the inner disk velocities, with the modeled pv-diagrams for two different velocities fields. On the left panel, we consider an infalling disk around a 20 M star (colors), where the gas is moving radially at 70% of the corresponding free-fall velocity (v ff = 2GM /R 1/2 ). The disk extends from the dust sublimation radius, where the disk temperature approaches 1500 K, up to a radius of 3000 au from the central star. We find that, this fraction of the free-fall velocity, combined with an inner cutoff near 500 au, best match the velocity range covered by the observed pv-diagram (green contours). Note that, since the infalling motion starts with non-zero velocities at the outer radius, this produces the inner hole near zero offsets in the modeled pv-diagram. The infalling profile well reproduces the observed pv-diagram in the second and fourth quadrants, which are the quadrants forbidden under the assumption of purely rotational motion (right panel).
In the right panel of Fig. 6, we show the comparison of the observed pv-diagram (green contours) with that expected for a disk in Keplerian rotation around a 20 M star (colors). We also draw the outer contour of the infalling disk model (cyan contour), in order to highlight the regions excluded by simple rotation. At variance with the infalling profile, we find that a Keplerian-like profile better reproduces the higher velocities in the first and third quadrants, when we assume an inner cutoff near 250 au and outer disk radius of 2400 au, similar to the deconvolved size of our disk profile (Fig. 1).
This analysis shows that the velocity field of gas through the disk is a combination of infalling and rotational motions, within a radius of 3000 au from the central star. On the one hand, the detection of a radial flow towards the central star implies that the condition of centrifugal equilibrium does not hold outside a radius of 500 au, where the velocity field has to be a combination of sub-Keplerian rotation and infalling motion. This scenario resembles that of the sub-Keplerian infalling disks modeled by Seifried et al. (2011) under the presence of strong magnetic fields. Our model predicts a mass infall rate at a radius of 500 au of 6 × 10 −4 M yr −1 . On the other hand, our data suggest that the inward gas flow slows down in the inner disk regions (≤ 500 au), where the pv-diagram resembles that of a centrifugally supported disk. However, these inner regions are only sampled by two beams at the current resolution. We explicitly note that, although we fixed the star mass in our model, the sim-ple agreement between observed and modeled pv-diagrams provides indirect confirmation of the central mass determined from the bolometric luminosity.
Moreover, in the right panel of Fig. 6 it is evident that there is an excess of blueshifted emission near zero offsets. This emission coincides with the molecular jet detected in our channel maps (right column of Fig. 4). If the molecular species does not trace the disk gas exclusively, and the disk plane is not seen edgeon, we show that the blueshifted outflow emission close to the star leaves a strong imprint in the pv-diagrams. In the Fig. 7, we provide the same pv-analysis for the CH 3 OH (18 3,16 -17 4,13 ) A + line emission. Interestingly, at the higher excitation-energy of this line, the blueshifted jet emission strongly affects the eastern side of the pv-diagram. In the channel maps (right panels), the spatial distribution of the jet emission bends to the redshifted side of the disk indeed, and progressively aligns to the outflow axis moving away from the disk plane. These maps definitely underline that the jet emission arises in the inner disk regions (< 800 au), providing a mechanism to transfer angular momentum away from the disk, and ensuring an inward flow of mass towards the central star.
Under steady state accretion, the mass infall rate of our model would imply that the final star mass could be as large as three times the current value, assuming the accretion phase lasts for about 10 5 yr. Nevertheless, it has been recently shown that the accretion process of young massive stars undergoes episodic accretion bursts, with the accretion rate suddenly rising by a few orders of magnitude (Caratti o Garatti et al. 2017;Hunter et al. 2017). Therefore, the net mass accretion onto the central star might exceed that inferred from the model.

Disk mass
In Fig. 8, we overlap a map of the observed dust continuum emission at 1.37 mm (dashed red contours) with that of the modeled dust emission from the disk at the same wavelength (white contours). We draw the contours of the observed continuum emission down to a 20% level of the peak of 14.6 mJy beam −1 . Below this level, corresponding to a radius of about 2000 au from the star, the dust continuum emission suffers contamination from the outer envelope, which is outlined by the lower, dashed, white contour in Fig. 1 (corresponding to about the 10% peak level). Our model predicts a dust continuum peak of 12.4 mJy beam −1 , in excellent agreement with the observed intensity.
We compare the disk mass estimated from the dust continuum emission with that predicted by the model. Following Hildebrand (1983), we estimate the disk mass from the continuum flux of 62.9 mJy within the 20% contour. We use an average gas temperature of 279 K, determined from the model between radii of 10 au and 2000 au, and assume a standard gas-to-dust mass ratio of 100 and a dust opacity of 1.0 cm 2 g −1 at 1.3 mm (Ossenkopf & Henning 1994, same for the model). From these parameters, we calculate a disk mass of 1.6 M and assign an uncertainty of ± 0.3 M , corresponding to an uncertainty of ± 20% of the measured flux. Our model predicts a disk mass of 2.5 M within a radius of 3000 au, which decreases to 1.4 M at 2000 au. Although the model does not account for a dusty envelope surrounding the disk, which might contribute to the observed dust flux, these values are in excellent agreement with the measured disk mass. We conclude that very likely the disk mass is much less (∼ 10%) than the mass of the star. Johnston et al. (2015) Previous ALMA observations by Johnston et al. (2015) reported about a circumstellar disk around a young star in AFGL 4176, which has a mass comparable to G023.01−00.41. These two sources have similar distances from the Sun and were observed with comparable resolution as well, allowing for a direct comparison. They might represent two snapshots for the formation of an O-type star, with G023.01−00.41 being younger, and having the potential to become three-times more massive than AFGL 4176 at the end of the accretion phase. In the following, we highlight common features and differences:

Comparison with
-Both disks extend up to radii of ∼ 2000 au from the central star, and have similar average temperatures of ∼ 200 K, as determined by LTE analysis of the CH 3 CN K-ladders. -There is a large difference (5×) in the estimate of the disk mass for these two objects, of 1.6 M and 8 M for G023.01−00.41 and AFGL 4176, respectively. However, this difference is due to the different dust opacities at 1 mm assumed in the calculations, of 1.0 cm 2 g −1 and 0.24 cm 2 g −1 respectively, whereas the integrated fluxes at the same wavelength are similar (50-60 mJy). We argue that the lower dust opacity used by Johnston et al. is more appropriate for diffuse clouds (Draine 2003). If consistent dust opacities are used, the two disks also have similar masses of ∼ 2 M , which amount to ∼ 10% of the central star mass. -The main difference between the two disks is related to the disk kinematics. On the one hand, Johnston et al. found Keplerian-like rotation up to the outer disk radii, meaning that nearly centrifugal equilibrium slows down any inward flow of mass through the disk. On the other hand, we find that the disk surrounding G023.01−00.41 is infalling and sub-Keplerian, and only at radii of a few 100 au might approach centrifugal equilibrium. According to Kuiper et al. (2011, their Fig. 4), larger Keplerian disks are expected at later times, whereas Keplerian rotation around the youngest stars is confined to the inner disk regions. This evidence suggests that G023.01−00.41 is still in an active phase of accretion and much younger than AFGL 4176.

Conclusions
We report about Atacama Large Millimeter/submillimeter Array (ALMA) observations, at wavelengths near 1 mm, of dense molecular gas and dust in the vicinity of an O-type young star, with a linear resolution as good as 700 au and a line sensitivity of 1 K. We targeted the luminous hot-molecular-core G023.01−00.41 after selecting the best candidate disk-jet system from our previous observations. We have resolved a (molecular) disk-jet system around a young star which currently attains a mass of 20 M . We present a kinematic analysis of the position-velocity diagrams of dense gas along the disk midplane, and compare them with the position-velocity diagrams simulated with a radiation transfer model. We show that the disk is falling in close to free-fall and slowly rotating with sub-Keplerian velocities, from radii of about 2000 au and down to 500 au from the central star, where we estimate a mass infall rate of 6 × 10 −4 M yr −1 . The disk mass is a small fraction of the star mass (∼ 10%). Furthermore, we are able to image the jet emission which arises from the inner disk radii (< 800 au), and show that its blueshifted emission leaves a strong imprint in the position-velocity diagrams.  Fig. 4, but for the CH 3 OH (18 3,16 -17 4,13 ) A + line emission with E up of 446 K. The right panels show two channel maps of the jet emission at a V LSR of +73.2 (middle) and +72.5 km s −1 (right). Fig. 8. Analysis of the dust continuum emission from the disk. Comparison of the dust emission observed with ALMA towards G023.01−00.41 (red dashed contours) with the modeled dust emission from a circumstellar disk around a 20 M star (colors and white contours). The brightness scale of the modeled emission is quantified by the righthand wedge. The red dashed contours are plotted at steps of 10% of the observed peak emission, starting from 90%; the white contours draw the same absolute levels for the model. The reference system and symbols are the same used in the previous figures. The synthesized ALMA beams, for the observed (red filled circle) and modeled (white circle) dust continuum maps, are shown in the bottom left corner. ing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. A.S. gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) Priority Program 1573. RK and AK acknowledge funding from the Emmy Noether research group on "Accretion Flows and Feedback in Realistic Models of Massive Star Formation" funded by the German Research Foundation under grant no. KU 2849/3-1.