Digging into the Interior of Hot Cores with ALMA: Spiral Accretion into the High-mass Protostellar Core G336.01–0.82

We observed the high-mass star-forming core G336.01–0.82 at 1.3 mm and 0.″05 (∼150 au) angular resolution with the Atacama Large Millimeter/submillimeter Array (ALMA) as part of the Digging into the Interior of Hot Cores with ALMA survey. These high-resolution observations reveal two spiral streamers feeding a circumstellar disk at opposite sides in great detail. Molecular line emission from CH3OH shows velocity gradients along the streamers consistent with infall. Similarly, a flattened envelope model with rotation and infall implies a mass larger than 10 M ☉ for the central source and a centrifugal barrier of 300 au. The location of the centrifugal barrier is consistent with local peaks in the continuum emission. We argue that gas brought by the spiral streamers is accumulating at the centrifugal barrier, which can result in future accretion burst events. A total high infall rate of ∼4 × 10−4 M ☉ yr−1 is derived by matching models to the observed velocity gradient along the streamers. Their contribution accounts for 20%–50% the global infall rate of the core, indicating streamers play an important role in the formation of high-mass stars.


INTRODUCTION
Non-spherical accretion flows are predicted to be key for the formation of high-mass stars as these allow to overcome the radiation pressure from the star.Accretion flows can form due to instabilities in the envelope/disks of hot cores forming high-mass stars (Oliva Fernández-López et al. 2023).In cases like AFGL 4176 mm1, these features appear within partially unstable Keplerian-like disk (Johnston et al. 2020).While in other cases (e.g., W33A and IRAS 18089-1732;Maud et al. 2017;Izquierdo et al. 2018;Sanhueza et al. 2021), these are feeding circumstellar disks.Recently, the Digging into the Interior of Hot Cores with ALMA (DI-HCA) project, showed a tentative example feeding a binary system (G335.579-0.272MM1 ALMA1; Olguin et al. 2022).
In this letter, we study a system of two accretion flows feeding a disk embedded in the high-mass star-forming region G336.01-0.82(distance: 3.1 kpc; Urquhart et al. 2018) as revealed by DIHCA observations.These observations reveal in great detail the motion along the accretion flows, allowing for future comparisons to numerical simulations.

OBSERVATIONS
ALMA band 6 (220 GHz) observations were performed by the 12 m array in two array configurations (similar to C43-8 and C43-5) as part of DIHCA.The long baseline observations (91 -8500 m) were made during July 2019 with 45 antennas while the short baseline observations (15 -1300 m) were performed during November 2018 with 44 antennas.The observations followed the same spectral setup as in Olguin et al. (2021Olguin et al. ( , 2022)), that is four spectral windows of ∼1.8 GHz with a spectral resolution of ∼976 kHz (∼1.3 km s −1 ).The maximum recoverable scale (MRS) of the long and short baseline data are 0. ′′ 94 (∼2900 au) and 3. ′′ 9 (∼12000 au), respectively.The data were calibrated using CASA versions 5.4.0-70 and 5.6.1-8reduction pipelines (CASA Team et al. 2022) for long and short baseline observations, respectively.Phase and amplitude self-calibration was then performed with incremental steps of shorter time intervals.The data from both configurations were independently self-calibrated (see Appendix A for a discussion of this method), with shortest time steps for phase self-calibration of 20 and 15 s for extended and compact configurations, respectively.A single time step of 30 s was used for amplitude self-calibration.Clean parameters used for self-calibration are the same as those used for continuum imaging (see below).Amplitude selfcalibration solutions were applied only to the continuum data.In order to maximize the resolution with the chosen robust parameter, we imaged the extended baseline data alone, while in some specific cases we imaged the combined data to recover extended emission as described below.
The line-free continuum and continuum-subtracted long baseline visibilities were obtained using the proce-dure described in Olguin et al. (2021).The continuum was then imaged using the tclean task of CASA with the Hogbom deconvolver and Briggs weighting with a robust parameter of 0.5.Figure 1(a) shows the continuum image of the region.We achieved an angular resolution of 0. ′′ 063 × 0. ′′ 044 (195 au × 136 au; position angle PA=-3.3 • ) and a noise level of 55 µJy beam −1 .In addition, we produced a continuum map from the combined data with the same parameters above.This map has an angular resolution of 0. ′′ 08 × 0. ′′ 59 (PA=-31.4• ) and a noise level of 54 µJy beam −1 .
In order to resolve the kinematics we imaged the long baseline data for CH 3 OH J Ka,Kc = 18 3,15 − 17 4,14 A, v t = 0 (233.795666GHz, E u = 447 K).Similarly the long baseline data of HC 3 N J = 24 − 23 l = 1f, v 7 = 1 (219.1737567GHz, E u = 452 K) was imaged to map the molecular gas distribution around the central source where CH 3 OH becomes optically thick.On the other hand, the combined data of SiO J = 5 − 4 and CH 3 CN J = 12 − 11 K-ladder was imaged to recover extended emission for the study of outflows and avoid missing flux to estimate gas temperatures.The imaging was performed using the auto-masking tool YCLEAN (Contreras et al. 2018).Noise levels range between 1.9-2.6 mJy beam −1 per channel, with maximum residuals on or below the 2σ level.With exception of SiO, the imaging was done with the same parameters of the continuum.For SiO the multiscale deconvolver was used.

RESULTS
The continuum data in Figure 1(a) shows three continuum structures (local peaks) resolved by the long baseline observations.Their peaks are located roughly in a line with a PA of 125 • , with the east and west peaks separated by roughly 400 au from ALMA1.The brightest structure is ALMA1-WEST (6.5 mJy beam −1 ), while ALMA1 (5.0 mJy beam −1 ) and ALMA1-EAST (3.4 mJy beam −1 ) are slightly fainter.However, Figure 1(b) shows that HC 3 N emission is single peaked, with a peak position slightly north of ALMA1.ALMA1 also seems to be the source powering the outflow seen in SiO (Figure 1(c)).Therefore, we conclude that ALMA1-EAST and ALMA1-WEST are not protostars and their nature is be discussed below.The direction of the outflow (PA ≈ 35 and 205 • for the red and blue shifted lobes, respectively) is relatively perpendicular to the line crossing the three structures.Close to the base of the outflow, the SiO emission widens which may point to a widening of the outflow cavity or indicative of accretion shocks.A spiral-like feature is connected to ALMA1-WEST coming from the north, while a similar feature is less evident to the south of ALMA1-EAST.A contin-   uum lane is observed toward the east oriented north to south which seems to be in part the result of the outflow interacting with the larger scale gas reservoir (see also Appendix B).We did not find any kinematic evidence linking the continuum lane to the east with the northern spiral-like feature.
To determine the nature of these continuum structures we explored the emission from CH 3 OH.Figure 2 shows moment maps of the CH 3 OH J Ka,Kc = 18 3,15 − 17 4,14 A, v t = 0 transition, while example spectra can be found in Appendix B. The moment maps are calculated over a velocity range of ∼12.5 km s −1 centered on the line frequency.For the zeroth moment map we used a variable integration window of twice the line FWHM for data over 5σ rms with σ rms = 2.6 mJy beam −1 and the whole velocity range for the remaining pixels, while for first and second moment maps we consider only emission over 5σ rms .To determine the line central velocity, we use a systemic velocity of -47.2 km s −1 (determined from 13 CH 3 CN by Taniguchi et al. 2023).As shown in Figure 2(b) this velocity is consistent with the velocity at the position of ALMA1.This would indicate that ALMA1-WEST and ALMA1-EAST correspond to gas rotating around the central source.Indeed, Figure 3(a) shows the position-velocity (pv) map following the arrow in Figure 2(a) with a slit width of 0. ′′ 05 (5 pixels).This is consistent with a combination of rotation and infalling motions (see §4.1).
In addition to rotation, Figure 2(b) shows more clearly the extent of the spiral feature connected to ALMA1-EAST, and reveals a velocity gradient along the spiral features.Figures 2(b) and 3(b) show that the northern spiral is increasingly blue-shifted as it join ALMA1-WEST, while Figures 2(b) and 3(c) show that the southern spiral becomes increasingly red-shifted as it joins ALMA1-EAST.We estimate a velocity gradient by fitting a line to the intensity weighted average velocity at each offset position.Note that both gradients, and in particular the northern one, are not completely linear and the gradient is steeper closer to the source, which may be caused by the change of the projection angle along the streamer combined with a change in velocity as the gas approaches the central source.These gradients point toward acceleration of the inflow, which has also been observed in other sources (e.g., Beltrán et al. 2018;Sanhueza et al. 2021).Figure 3 (b) also shows a main spine associated to the spiral and fainter emission towards higher velocities, e.g., between offsets 0. ′′ 2 and 0. ′′ 6 over the systemic velocity dashed green line.This fainter emission may be the result of multiple velocity components caused by the envelope gas interacting with the spiral as it rotates or gas associated with the redshifted outflow lobe.Linewidth (km s 1 )  To interpret the CH 3 OH emission we use a infalling and rotating envelope (IRE) model as implemented in FERIA (Oya et al. 2022).The IRE model consists of a rotation velocity: and an infall velocity: with G the gravitational constant, M the central mass and r cb the radius of the centrifugal barrier.FERIA can also include a Keplerian disk, however the CH 3 OH transition in Figure 3(a) is more optically thick towards ALMA1, hence we fix the inner radius of the model r in = r cb .The IRE model relative intensity is determined from a combination of the contribution of the density profile (∝ r −1.5 for free-fall, e.g., Shu 1977) and the temperature profile (∝ r −0.4 for high-mass star forming regions, e.g., Gieser et al. 2021).FERIA convolves the model results with a Gaussian beam and produces pv maps.We computed a grid of models with inclination angles between 50-80  • and a radius of 800 au.Similar good fits are obtained for models with inclination angles between 60 and 75 • and r cb of 200 and 300 au.In general, a larger r cb favors a larger inclination angle, e.g., the best model with r cb = 300 au has an inclination angle of 75 • .These trends are relatively independent of the model mass and outer radius.Due to the symmetry of the model, differences in intensity between the ALMA1-EAST and ALMA2-WEST cannot be reproduced.Additionally, since FERIA works under the optically thin approximation, line profiles cannot be well reproduced towards the central region.Note that the PA of the line peak intensity of the model shown in Figure 4(a) does not exactly match the observed one.However, the PA from the peaks of the model zeroth moment map is close to the observed one (124 • vs. 125 • ).We expect a closer to edge-on model to have PA of the peak line closer to the one derived from the zeroth moment.Additionally, these spots may not be located at exact opposites as the model predicts due to the model's symmetry.Roughly, the velocity at the peak intensity in the pv maps constrains the central mass, while the offset of the peak constrains the centrifugal barrier (e.g., Oya et al. 2014).Similar to Csengeri et al. (2018), the peak positions of CH 3 OH are consistent with the radius of the centrifugal barrier (Figure 2(a)).
The position of the centrifugal barrier agrees with the continuum peaks ALMA1-EAST and ALMA1-WEST.As the gas falls into the disk through the spiral features, accretion shocks at the centrifugal barrier (e.g., in the hot corino IRAS 16293-2422 Source A; Oya et al. 2016Oya et al. , 2018;;Maureira et al. 2022) can result in accumulation of gas and enhancement of the gas temperature (e.g., Oya & Yamamoto 2020).Accumulation of matter at the location of the centrifugal barrier is predicted by simulations (Oliva & Kuiper 2020).Therefore, the continuum and line data seems to be showing an enhancement of such features.Without further high-resolution observation we can only speculate whether there are additional substructures connecting these features to the inner regions of the disk.

Streamer properties
The CH 3 OH lines with low E u tend to be optically thick precluding a fit for temperature determination, even toward the spiral arms.In order to estimate the gas temperature, we thus use the average brightness temperature of CH 3 CN J = 12 − 11 K = 2 to 4 transitions at peak emission.These CH 3 CN K-level transitions are chosen because they appear optically thick (line peaks converge toward a constant temperature) and are not blended with other lines.As lines become optically thick their brightness temperature approaches the gas temperature.The brightness temperature of the CH 3 CN are roughly the same between the combined and extended configuration, which implies that the emission fills the beam (the emission is extended), but also that not much is resolved out in the extended baseline data alone.Figure 5 shows the temperature map while Appendix C shows the standard deviation of the peak values.
The mass in the spirals is calculated from the continuum using: with S ν the continuum flux density, d the distance to the source, R dg = 0.01 the dust-to-gas ratio, κ 1.33mm = 1 cm 2 g −1 the dust opacity (Ossenkopf & Henning 1994), and ral masses, M sp , of 1.3 and 0.5 M ⊙ for the northern and southern spirals, respectively.

Infall along streamers
To study infall motions along the streamers, we computed models of gas collapsing along parabolic streamlines.The velocity distribution follows that of a rotating and collapsing envelope under angular momentum conservation (Mendoza et al. 2009), and is described by the central mass M , centrifugal radius r c , outer radius R out , inner radius R in , initial polar angle θ 0 and initial radial velocity v r0 .The streamlines are calculated taking into consideration the projection angles on the sky and the initial azimuthal angle ϕ 0 .The models were computed following the implementation of the equations by Pineda et al. (2020).We used a PA=125 • (see §3), and following the IRE model, we set M = 10 M ⊙ and R in = r c = 2r cb = 400 au. Figure 4 show that the streamer distribution and line of sight velocity are well described by the streamline models.The models have R out = 2500 au, while the initial values are θ 0 = 80 • , ϕ 0 = 60 • and 280 • , and v 0 = 0 and 2 km s −1 for the northern and southern spirals, respectively.
In order to estimate the average infall rate along the streamers, first we estimate the volume weighted average velocity from the discrete points evaluated by the model in Figure 4(b) as: ΣΣv(r, θ)r 2 sin θ∆r ∆θ ΣΣr 2 sin θ∆r ∆θ (4) where (r, θ) are the radial and polar angle in spherical coordinates.We obtain average infall velocities of 2.7 and 4.0 km s −1 for the northern and southern spirals, respectively, and a maximum velocity of ∼6 km s −1 at r c .Assuming a cylindrical streamer, the infall rate can be calculated as Ṁ = v inf M sp /l where l is the length of the streamer (e.g., Kirk et al. 2013).Using the masses from §4.2 and lengths of 2870 and 2370 au, we obtain infall rates of 2.5 and 1.8 × 10 −4 M ⊙ yr −1 for the northern and southern streamers, respectively.Alternatively, the infall rate can be calculated from the total time it takes a parcel of gas to go from r 0 to r c for each streamer, and also from the free-fall time.
For the former, we use the streamer trajectories and velocities to estimate infalling times of 7.6 kyr for the northern and 3.5 kyr for the southern streamers, yielding infall rates of 1.7 and 1.4 × 10 −4 M ⊙ yr −1 .On the other hand, we obtain a free-fall time, with M tot the total mass (streamers and central source), of 5.8 kyr.This in turn gives infall rates of 2.2 and 0.8 × 10 −4 M ⊙ yr −1 .These results indicate a contribution of roughly 1 − 2.5 × 10 −4 M ⊙ yr −1 from each streamer to the total accretion rate.
On the other hand, the core flux density is 36 mJy as measured in the compact configuration data (∼0.′′ 3 = 930 au resolution, MRS > 10 4 au) from dendrogram identification of the cores (Ishihara et al., in preparation).From eq. 3 we obtain core masses of M c = 1.1 − 2.4 M ⊙ assuming a dust temperature of 100 and 50 K, respectively.For a core of radius R = 2500 au (i.e., the same radius as the origin of the streamers), we estimate an average infall velocity (eq.4) v inf = 3 km s −1 based on the infall velocity distribution given by eq. 2. Hence for a spherically symmetric core the infall rate is Ṁ = 3v inf M c /R = 0.8 − 2 × 10 −3 M ⊙ yr −1 .This imply that the contribution of the streamers may account for a fifth to a half of the core infall rate.Therefore, their contribution plays a significant role in increasing the central source mass.
Infall streamers have been reported in low-mass (e.g., Pineda et al. 2020;Kido et al. 2023; see also Pineda et al. 2023 for a review) as well as in high-mass star formation.Their presence have been invoked to explain accretion bursts, resulting in infall/accretion rates and luminosity higher than expected (e.g., Valdivia-Mena et al. 2022).This is also supported by simulations, which in the case of high-mass star-forming regions predict bursts of high-accretion rates as the gas from the streamer fragments is consumed (e.g., Meyer et al. 2018).Our observations support this scenario, even though multi-wavelength and/or multi-epoch observations are still needed to trace changes in luminosity.However, streamers in isolated star formation simulations seem to dominate scales below 1000 au (e.g., Meyer et al. 2018;Oliva & Kuiper 2020), while observations like those presented here have shown that these features can appear well beyond the disk radius (defined as the centrifugal barrier).As isolated cores, these simulations do not accurately reflect the relation between inner regions of the cores and inflow(s) from the larger clump scales.On the other hand, simulations considering interactions between cores can explain observations of asymmetric streamers at > 1000 au scales (e.g., Kinoshita & Nakamura 2022, Yano et al., submitted), like in W33A where a single spiral feature is observed.Therefore, whether the streamers in G336.01-0.82originated from gravitational instabilities, large scale inflows and/or interaction with other cores remain unclear.The lifetimes of these structures can only be answered by large samples like the one provided by the DIHCA survey.

CONCLUSIONS
We present ALMA observations of the high-mass core G336.01-0.82 at 1.3 mm with an angular resolution of 0. ′′ 05 (155 au).The continuum data shows three peaks with two of them connected to spiral-like structures or streamers.Velocity gradients are observed along the streamers in CH 3 OH line emission, suggestive of infall, while SiO emission shows a bipolar outflow originating from the central source.An infalling and rotating model indicates a central source mass of 10 M ⊙ and a centrifugal barrier of 300 au, consistent with the position of the continuum peaks.We estimate the gas temperature from the brightness temperature of four CH 3 CN transitions.From the temperature map and the continuum emission we derive masses of 1.3 and 0.5 M ⊙ for the northern and southern spirals, which results in a total accretion rate of ∼4 × 10 −4 M ⊙ yr −1 based on streamline models.The position of the dust continuum peaks seem to indicate that the gas being fed by the streamers is accumulating at the outskirts of the disk.This may originate future accretion burst events.The streamers also could be unstable which would lead to fragmentation and eventually the formation of stellar companions.
Given that DIHCA data was delivered during different cycles for each array configuration, self-calibration was performed on individual configurations in order to optimize time.However, a better approach would be to perform self-calibration on the combined data sets.As such, we tested the changes on images produced by these approaches with phase only self-calibration tables.For these tests we produced continuum images in a uniform way: continuum visibilities (with the same line free channels) were CLEAN with the same fixed mask and thresholds.In order to account for shifts of the baselines of the individually self-calibrated images, we also performed an additional self-calibration correction.This correction was obtained by using the model of the last cleaning step of the extended configuration self-calibration and an infinite solution interval.We compared extended configuration alone and combined data.Overall, we obtained changes in intensity and flux density (over the same central area) of < 5%, with pixels with lower intensities accounting for the highest changes.The results for extended configuration alone are independent of whether the additional self-calibration step to align the baselines is performed.
Additionally, we calculated dirty cubes of a subset of channels where lines are expected.We found changes in intensity of < 2% and no significant changes in line profiles.In summary, the two approaches account for less than the typical ALMA flux calibration error of 10%, and part of the changes can also be attributed to the different CLEAN steps performed during self-calibration.However, we would like to note the experiment described above was performed on a typical DIHCA target with high S/N ratio (≳100).For sources with lower signal-to-noise ratios differences may be obtained, as combining the data sets would increase the S/N ratio allowing to make better models for self-calibration.C. TEMPERATURE ERROR MAP Figure 8 shows the standard deviation of the peak brightness temperature of CH 3 CN J = 12 − 11 K = 2 to 4.

Figure 2 .Figure 3 .
Figure 2. CH3OH JK a ,Kc = 183,15 − 174,14 A, vt = 0 moment maps.(a) shows the zeroth order moment map.The blue and red lines mark the position of the north and south spiral-like features, respectively.The arrow marks the direction of the pv map in Figure 3(a), and corresponds to the line crossing all three continuum structures (PA=125 • ).(b) shows the first order moment map, and (c) the second order moment map.The continuum is shown in contours with the same levels as in Figure 1, and the continuum peaks are marked with green triangles.The beam size ellipse is shown in the lower left corner.

Figure 4 .
Figure 4. Models explaining the kinematics of gas around G336.01-0.82from CH3OH JK a ,Kc = 183,15 − 174,14 A, vt = 0 emission in position-position-velocity space.This figure is available online as an interactive figure.The blue surface layers show the CH3OH emission levels of 5 − 15 × σrms in steps of 2.5 with σrms = 2.6 mJy beam −1 .(a) Shows the IRE model surface layers at the same emission levels.(b) Shows the spatial projection of the streamer models.In the interactive version, the velocity dimension can be explored and the models plotted enabled/disabled.The coordinates are with respect to ALMA1 (indicated with a yellow sphere: 16 h 35 m 09.s 261, -48 • 46 ′ 47. ′′ 659, -47.2 km s −1 ).The green lines indicate the projection of the zero offsets on each axis.
(5 • steps), central masses of 5, 8, 10, 15 and 20 M ⊙ , and r cb between 100-600 au (100 au steps) for IREs with 700, 800, 900 and 1000 au radii to match the extension of the CH 3 OH emission.For each model, χ 2 values were calculated from the pv map data over 3σ.The best-fitting model pv map is shown in Figure 3(a) and its 3-D distribution is shown in Figure 4.This model has a mass M = 10 M ⊙ , r cb = 200 au, inclination of 65 Figure 5. Average peak brightness temperature from CH3CN J = 12 − 11 K transitions 2 to 4. The blue and red regions in (a) indicate where the properties in §4.2 are estimated, and roughly follow the continuum 5σ level.Gray contours show the continuum emission from the combined data at the same levels as in Figure 1 but with σrms = 54 µJy beam −1 and the green triangles indicate the continuum peaks.The beam size ellipse is shown in the lower left corner.