Disk Kinematics and Stability in High-Mass Star Formation: Linking Simulations and Observations

In the disk-mediated accretion scenario for the formation of the most massive stars, gravitational instabilities in the disk can force it to fragment. We investigate the effects of inclination and spatial resolution on observable kinematics and stability of disks in high-mass star formation. We study a high-resolution 3D radiation-hydrodynamic simulation that leads to the fragmentation of a massive disk. Using RADMC-3D we produce 1.3 mm continuum and CH3CN line cubes at different inclinations. The model is set to different distances and synthetic observations are created for ALMA at ~80 mas resolution and NOEMA at ~0.3''. The synthetic ALMA observations resolve all fragments and their kinematics well. The synthetic NOEMA observations at 800 pc (~300 au resolution) are able to resolve the fragments, while at 2000 pc (~800 au resolution) only a single slightly elongated structure is observed. The position-velocity (PV) plots show the differential rotation of material best in the edge-on views. As the observations become less resolved, the inner high-velocity components of the disk become blended with the envelope and the PV plots resemble rigid-body-like rotation. Protostellar mass estimates from PV plots of poorly resolved observations are therefore overestimated. We fit the emission of CH3CN lines and produce maps of gas temperature with values in the range of 100-300 K. Studying the Toomre stability of the disks in the resolved observations, we find Q values below the critical value for stability against gravitational collapse at the positions of the fragments and the arms connecting the fragments. For the poorly resolved observations we find low Q values in the outskirts of the disk. Therefore we are able to predict that the disk is unstable and fragmenting even in poorly resolved observations. This conclusion is true regardless of knowledge about the inclination of the disk.


Introduction
High-mass stars (M 8 M ) live short and violent lives, ejecting a significant amount of mechanical and radiative energy into the Universe, enriching the interstellar medium with the heavy material needed for the creation of organic life as we know it. They affect not only subsequent star and planet formation, but also the evolution of the galaxy. Despite their importance, the short-lived nature of their existence makes studying their birth and evolution extremely difficult. To make matters worse, high-mass stars are rare (Kroupa 2001;Chabrier 2003) and typically found at distances much further than low-mass star-forming regions.
High-mass stars experience rapid Kelvin-Helmholtz contractions, allowing them to start burning hydrogen while still accreting material (Palla & Stahler 1993). This lack of pre-mainsequence phase in the evolution of high-mass stars would mean that intense (ionising) radiation from the protostar is present in the main accretion phase of the protostar. For many years, the effect of radiation pressure on dust was thought to hinder the formation of stars more massive than ∼40 M (e.g., Larson & Starrfield 1971;Kahn 1974;Wolfire & Cassinelli 1987). Such studies would adopt nearly dust-free cloud conditions and extremely high accretion rates to allow the build-up of enough material onto the protostar to even reach such stellar upper mass limits. Subsequent works showed that the unrealistic restrictions on the dust opacity and accretion rates can be removed if nonspherical accretion flows are considered (e.g., Nakano 1989). In particular, magnetic fields, rotation or even simple contraction of the cloud could provide such deviations from spherical symmetry. In more recent years, two-and three-dimensional radiation-hydrodynamical simulations of collapsing cores have shown support for disk-mediated accretion in the formation of very massive stars, analogous to the formation of low-mass stars (e.g., Yorke & Sonnhalter 2002;Krumholz et al. 2009;Kuiper et al. , 2011Kuiper & Yorke 2013;Klassen et al. 2016;Rosen et al. 2016;Kuiper & Hosokawa 2018). In the diskmediated accretion, the effect of radiation pressure would be re-Article number, page 1 of 20 arXiv:1909.04051v2 [astro-ph.GA] 15 Oct 2019 A&A proofs: manuscript no. sims_arXiv_v2 duced as the radiation could escape through the poles along the disk rotation axis, and the disk would be shielded from radiation due to its high optical depth.
While the formation mechanism of high-mass stars through accretion disks is being comprehensively studied theoretically, the observational existence of disks around O-and B-type stars is still elusive. This is in part due to their scarcity, fast evolution, and the fact that high-mass stars form in dense distant clusters, making disentangling individual stellar contributions difficult (see review by Motte et al. 2018). Perhaps the best observational evidence for the existence of disks is the ubiquitous observations of collimated bipolar outflows (e.g., Beuther et al. 2002;Fallscheer et al. 2009;Leurini et al. 2011;Frank et al. 2014;Maud et al. 2015;Sanna et al. 2019), which has also been predicted by theoretical models (e.g., Pudritz et al. 2007;Banerjee & Pudritz 2008;Kölligan & Kuiper 2018). A review by Beltrán & de Wit (2016) summarises our current understanding of the properties of accretion disks around young intermediate-to high-mass stars. They refer to large structures for which a rotational velocity field could be established, extending thousands of au in space, as 'toroids' which host a (proto)cluster of stars. Furthermore, with the advent of long-baseline interferometry, a hand-full of disk candidates around the most massive stars have recently been detected (Johnston et al. 2015;Zapata et al. 2015Zapata et al. , 2019Ilee et al. 2016Ilee et al. , 2018Chen et al. 2016;Cesaroni et al. 2017;Csengeri et al. 2018;Ahmadi et al. 2018;Maud et al. 2018Maud et al. , 2019Sanna et al. 2019).
Although disk-mediated accretion has solved the radiationpressure problem, high accretion rates ( 10 −4 M yr −1 ) through the disk are still required for the creation of the most massive stars. Gas densities needed to provide such high rates of accretion could induce gravitational instabilities in the disk forcing it to fragment and produce companion objects (Kratter & Matzner 2006). While high densities can induce instabilities in the disk, thermal gas pressure and the shear force as a result of differential rotation of material in the disk can provide added stability against local collapse. The balance between these forces ultimately determines the fate of the disk. Originally introduced by Safronov (1960) and later quantified by Toomre (1964), a disk in Keplerian rotation is unstable against axisymmetric gravitational instabilities when the Toomre Q parameter where the stabilising effect of pressure is accounted for in the equation of sound speed c s , shear is considered in the epicyclic frequency (or angular velocity) Ω of the disk, and the selfgravitational force is represented as the surface density Σ of the disk. Disks are prone to fragmentation when their mass exceeds 50% of the mass of the entire system (Kratter et al. 2010). Numerical simulations have found that disk fragmentation can result in the creation of short-period binaries on the scales of hundreds of au (Meyer et al. 2017;. The highest-resolution observations with the Atacama Large Millimeter Array (ALMA) that are able to resolve such structures on few hundred au scales are also starting to see fragmentation of the disks on these scales (Ilee et al. 2018). Considering the large fraction (> 70%) of OB stars that are found in close binary systems through radial velocity surveys (Chini et al. 2012), it is important to understand the role of disk fragmentation in affecting the observed stellar mass distribution.
As we obtain the observational capabilities to resolve highmass protostellar disks, it is important to create frameworks with which to study their properties. Synthetic observations play a crucial role in this effort as they not only allow for meaningful predictions to be made from theoretical models, but also provide valid interpretations for observations (for a review, see Haworth et al. 2018). To date, only a few works have created synthetic observations from radiation-hydrodynamic simulations of high-mass star formation with the aim of providing predictions for long-baseline interferometry (Krumholz et al. 2007;Harries et al. 2017;Meyer et al. 2018). Due to the computational cost of running numerical simulations and radiative transfer models for a wide range of parameters, such studies have not been very comprehensive. Jankovic et al. (2019) were able to explore a larger parameter space by adopting a semi-analytical treatment of the density, temperature, and velocity fields with prescriptions for molecular freeze-out and thermal dissociation. Their synthetic dust continuum and molecular line observations probe different disk masses, distances, inclinations, thermal structures, dust distributions, and number and orientation of spirals and fragments.
In the following work, we aim to study the kinematics and stability of high-mass protostellar disks in detail by creating synthetic observations for the highest resolution 3D radiationhydrodynamic simulations that lead to the fragmentation of a massive disk. In particular, we study the effect of inclination and spatial resolution on the derived disk and protostellar properties. To investigate the effects of spatial resolution, we assume a close distance of 800 pc and a more typical distance of 2 kpc, and create synthetic observations in the millimeter regime for ALMA at long baselines reaching very high resolutions (∼80 mas) and for the NOrthern Extended Millimeter Array (NOEMA) with 6 antennas covering much shorter baselines currently providing the highest resolutions possible (∼0.3 ) in the northern sky. In this way, we can study the properties of high-mass protostellar disks when they are resolved as well as unresolved disks which resemble toroidal-like structures often seen in observations of intermediate-to high-mass stars. The motivation for this work is in part because we have observed the largest sample of 20 highmass young stellar objects in the northern sky with NOEMA ) 1 and aim to study the stability of the candidate disks and characterise their properties in a follow-up paper. It is therefore paramount to understand the observational limitations in detail and benchmark our methods. From the analysis of synthetic observations, not only do we learn about the observational limitations but will also be able to use the observations to their full potential and gain higher confidence in the conclusions.
The paper is organised as follows. Sect. 2 outlines the details of the simulations. The radiative transfer modelling scheme is summarised in Sect. 3. In Sect. 4 we describe how the synthetic observations for ALMA and NOEMA were created. In Sect. 5 we summarise our analysis and results, divided into subsections describing both the continuum (Sect. 5.1) and gas kinematics (Sect. 5.2), temperature distribution (Sect. 5.3), mass estimates from dust and position-velocity diagrams (Sect. 5.4), and the Toomre stability of the disks (Sect. 5.5). A summary of the work and main conclusions can be found in Sect. 6.

Simulations
The numerical simulation starts from a collapsing reservoir of dusty gas and follows the evolution of the system including the formation and evolution of a central massive protostar and the formation and fragmentation of its surrounding accretion disk.
A theoretical investigation of this simulation (plus a convergence study) will be presented in Oliva et al. (in prep.) and covers an analysis of the disk's stability against non-axially symmetric features such as spiral arms, an analysis of its fragmentation probability, a determination of the fragment statistics, their masses, and orbital parameters. For the purpose of the synthetic observations presented herein, we summarise in the following the covered physics in the simulation, the limitations of the model, the physical initial conditions, and the numerical grid setup, highlighting the uniqueness of the extremely high spatial resolution of the disk fragmentation simulation. We refer to Kuiper et al. (2011), Meyer et al. (2017), and Meyer et al. (2018 for previous three-dimensional studies of disk evolution around highmass protostars that use the same numerical framework. For the hydrodynamic evolution, we utilise version 4.1. of the open source code Pluto (Mignone et al. 2007;Mignone et al. 2012). Self-gravity of the gas is included via the Haumea module Kuiper et al. 2011). Stellar irradiation and disk cooling are modelled using the hybrid radiation transport module Makemake Kuiper & Klessen 2013), updated to a so-called two-temperature scheme (Commerçon et al. 2011). The first application of this scheme is presented in Kuiper & Hosokawa (2018) and Bhandare et al. (2018). The evolution of the central protostar follows the tracks by Hosokawa & Omukai (2009). Details on the numerical implementation and coupling to the hydrodynamics are given in  and Kuiper & Yorke (2013).
The simulation focusses on the formation, evolution, and fragmentation of the circumstellar accretion disk. In order to be able to achieve the required high spatial resolution (discussed in detail below), we have to neglect other potential feedback effects from the forming protostar. First, radiation forces are not taken into account. Due to the fact that we analyse a snapshot of the simulation, in which the central protostar has reached 10 M , radiation forces are indeed much lower than gravity and centrifugal forces. Second, photoionisation feedback and the formation of an HII region is not taken into account. Again, those physics take place at a later evolutionary stage than considered here, namely after the contraction of the central protostar to the zero age main sequence, when it starts to emit copious radiation in the EUV regime. With respect to feedback of radiation and pressure and photoionisation, we refer to Kuiper & Hosokawa (2018). Third, no feedback by protostellar outflows is taken into account in this simulation unlike in our previous feedback studies (Kuiper et al. 2015;Kuiper et al. 2016). The injection of a fast outflow would reduce the numerical timestep per main iteration due to the CFL condition (Courant et al. 1967) of the explicitly solved hydrodynamics. This would make the simulation computationally too expensive and modelling of the entire disk fragmentation phase would become infeasible. From a physical point of view, the outflow itself will only marginally impact the disk physics by lowering the optical depth in the bipolar outflow cavity. The most severe limitation of the current disk fragmentation model is the negligence of large-scale magnetic fields pervading the initial collapsing mass reservoir. Although we have previously used the same numerical framework to model the collapse of magnetised high-mass star forming regions (Kölligan & Kuiper 2018), the effect of magnetic fields are neglected here. The collapse of a magnetised mass reservoir would not only lead to the launching of fast, collimated jets and wide-angle outflows, but actually also impact the morphology and dynamics of especially the innermost disk structure via additional magnetic pressure and angular momentum transport. The final impact of magnetic fields on disk fragmentation has to be investigated in future studies.
The initial mass reservoir of the simulation is described by a spherically symmetric mass density profile ρ(r) ∝ r −1.5 and an axially symmetric rotation profile Ω(R) ∝ R −0.75 , where r and R denote the spherical and cylindrical radii, respectively. The outer radius of the reservoir is set to 0.1 pc and its total mass to 200 M . The corresponding free-fall time is t ff ≈ 37 kyr. The ratio of rotational energy to gravitational energy of the initial mass reservoir is set to 5%, corresponding to an angular rotation frequency normalised to a value of roughly 10 −10 Hz at R = 10 au. The temperature of the gas and dust is initially set to 10 K.
The uniqueness of this simulation is based on its extremely high spatial resolution of the forming accretion disk and its fragments. The gravito-radiation-hydrodynamical evolution of the system is computed on a numerical grid in spherical coordinates. In order to achieve a high dynamic range, we refine the grid cells in the radial direction toward the inner disk rim and in the polar direction toward the disk midplane. The refinement is done such that the grid cells in the disk midplane (at θ = 90 • ) have the same extent in the radial, the polar, and the azimuthal direction (∆r = r∆θ = r∆φ). The radial extent of the inner sink cell is 30 au. The outer radius of the computational domain is 0.1 pc. In the polar direction, the grid extends from the upper to the lower polar axis (0 • ≤ θ ≤ 180 • ). In the azimuthal direction, the grid covers the full angle range as well (0 • ≤ φ ≤ 360 • ). The number of grid cells in the radial, the polar, and the azimuthal directions are 268 × 81 × 256, respectively. This grid configuration results in a spatial resolution of the disk's midplane layer as function of distance to the host star r of ∆x(r) ≈ 0.025 au × r/(1 au). (2) Hence, at the inner rim of the computational domain at 30 au, the disk is resolved down to 0.75 au. We analyse the disk structure at 12 kyr after the onset of collapse when the disk has grown to ≈ 1000 au (see Fig. 1). At this outer disk edge, the disk is resolved down to 25 au. In conclusion, the circumstellar disk (defined as gas with densities higher than 10 −16 g cm −3 ) is represented by more than a million grid cells in total. Both important physical length scales associated with disk evolution, namely the disk's pressure scale height H (which is related to disk cooling) and the disk's Jeans length λ J (which is related to the fragmentation process of the disk) are resolved. In comparison, to achieve the same dynamic range on a Cartesian grid with adaptive-meshrefinement or nested levels, 15 levels of refinement would be required. This is calculated as follows: to resolve the entire computational domain with L box = 2 × 0.1 pc down to a minimum spatial resolution of dx = 0.75 au requires 2 n = L box /dx, so that the levels of refinement n ∼ 15.7. Besides the high spatial resolution, another important and unique feature of this simulation is that no subgrid models such as sink particles are used to represent fragments in the disk. A strongly self-gravitating subregion of the disk can undergo local collapse, and fragmentation barriers such as the required disk cooling and disk shear timescale are self-consistently taken into account due to the fact that the model solves for the hydrodynamics, radiation transport, and self-gravity, while resolving the pressure scale height and Jeans length.The high spatial resolution of the grid allows us to even resolve the hydrostatic physics of the forming first Larson cores, at least in the inner parts of the disk where fragmentation occurs first; fragments, which are expelled from the disk toward larger radii (r > 1000 au) are not resolved first Larson cores anymore. Those first Larson cores, where self-gravity is in balance with internal thermal pressure, further interact with the large-scale accretion disk in which they are embedded. They accrete gas and dust from their own smallerscale accretion disks (discussed later in Sect. 5.2). The fragments migrate due to angular momentum exchange with the disk via gravitational torques, and sometimes merge, losing parts of their outer atmospheres back into the disk phase. Some fragments even get fully disrupted in the case of, for example, interactions with a spiral arm.
From this simulation, we take a snapshot at 12 kyr after the onset of collapse, when the host star has grown to 10 M . While at this snapshot the protostar would be classified as a B-type star, the setup of the simulation is valid for the formation of a more massive O-type star as it will continue to accrete through the disk and the migration of the fragments which induces episodic accretion onto the central protostar. Therefore, we expect this snapshot to also be representative of a slightly later evolutionary phase in which the protostar will be classified as an O-type star. In the simulations of Meyer et al. (2018) which follow a similar numerical setup as ours and reach stellar masses as high as 35 M , the effects of stellar irradiation and the initial angular velocity distribution are tested. In all cases, the disk forms spiral arms that fragment due to gravitational instabilities as in the simulation analysed in this work. What the change in the initial angular velocity distribution affects is the disk size and the temporal evolution of the disk to star mass ratio, but the fact that the disk would eventually fragment remains unchanged. Similarly, the adoption of a steeper density distribution as done in some other works (e.g., ρ(r) ∝ r −2 in Yorke 2013 andHosokawa 2018) would make the disk initially more stable as more mass would be concentrated in the center rather than distributed in the envelope/disk and this lower disk-to-star mass ratio would make the disk more Toomre-stable initially, but not enough to hinder spiral arm formation and ultimately the fragmentation of the disk (Kuiper et al. 2011). Therefore, the initial conditions we have chosen for this simulation are reasonable for the formation of massive protostar(s) and the study of disk fragmentation. Figure 1 shows the surface density map extracted from the numerical simulations and inclined to 10 • , 30 • , 60 • , and 80 • , where 0 • denotes a face-on view and 90 • denotes an edge-on view of the circumstellar disk. To extract the surface density map, we have made use of the pressure scale height where c s is the sound speed and v KSG the 'Keplerian' velocity taking into account the central host mass as well as the disk mass interior to that radius R. The surface density is then calculated via where ρ(x, y) is the gas mass density of the midplane. Since the calculation of the pressure scale height is only meaningful inside the disk region, the high column densities seen in the outskirts of the maps which are especially visible in the 80 • inclined map are not valid. The fragmented accretion disk of the snapshot at 12 kyr is further analysed by radiative transfer post-processing, as described in the following section.

Radiative Transfer Post-Processing
The simulation snapshot under investigation consists of a central protostar with 10 M , a radius of ≈ 9 R , and an effective surface temperature of 17 700 K, surrounded by a fragmented We will present the generated models in subsequent sections and in the following paragraph summarise the parameters used in RADMC-3D.
Gas and dust densities, temperature, and gas velocity are extracted from the hydrodynamics simulation snapshot. Gas and dust temperature are assumed to be in equilibrium, as it is done for the radiation-hydrodynamics simulation. Dust opacities are taken from Laor & Draine (1993), used also for the continuum radiation transport in the hydrodynamics simulation. The CH 3 CN abundance (with respect to H 2 ) is set to 10 −10 for temperatures ≤90 K, 5 × 10 −9 for temperatures >90 K, and 10 −8 for temperatures >100 K based on Collings et al. (2004) and Gerner et al. (2014). The spectral width of the ray-tracer is set to ±247.8 km s −1 around the central frequency of the non-Doppler-shifted line profile. The spectral resolution is set to 1400 frequency bins. The spatial resolution is 400 × 400 pixels for a region of the inner 8400 au × 8400 au around the protostar. The inclination-dependent output of this radiative transfer postprocessing task, mimicking an observation by a perfect telescope, is further processed to generate synthetic observational data cubes. Details are given in the following section.

Synthetic Observations
We shift the inclined post-processed snapshot of the model simulations at 12 kyr in both dust continuum and lines to two different distances of 800 and 2000 pc by dividing the fluxes by a factor of distance squared and adjusting the pixel resolution accordingly. The generation of the synthetic observations for ALMA and NOEMA are described in detail in the following subsections. The rms noise in the emission-free region. (b) The rms noise in the emission-free region in the channel that has the peak of emission for CH 3 CN (12 4 − 11 4 ), therefore the unit is per 0.35 km s −1 channel.

ALMA
To create synthetic ALMA observations, we use the ALMA simulator task simalma within the Common Astronomy Software Applications (CASA) software package 2 to create visibility files. The model simulations adjusted to the desired distance and in units of Jansky per pixel are used as the 'sky model'. The location of the sky model is set to a generic position in the southern sky, with the coordinates of the center of the map at α(J2000) = 24 h 00 m 00 s , δ(J2000) = -35 d 00 h 00 s . The observed sky frequency is set to 220.679 GHz with the widths of the channels for the line cubes set to 260.3 MHz (∼0.35 km s −1 ). We use a combination of two ALMA Cycle 5 configurations, 5-6 and 5-9, to reach a desired angular resolution of ∼0.08 . The configuration files containing the information about the distribution of the antennas were obtained within CASA (version 5.1.2-4) and were labeled 'alma.cycle5.6.cfg' and 'alma.cycle5.9.cfg' for configurations 5-6 and 5-9, respectively. In Cycle 5, ALMA made use of 43 antennas in the 12-m arrays. The uv-coverage of the ALMA array in configuration 5-6 has baselines in the range of 15 m to 2.5 km, and in configuration 5-9 covers a range of 375 m to 13.9 km. We set the integration time of the synthetic observations to 30-second intervals for a total time of 18 minutes in configuration 5-6 and 35 minutes in configuration 5-9. Figure 2 shows the uv-coverage of the synthetic observations. Thermal noise was automatically added to the observations assuming precipitable water vapour of 0.5 mm. The images were generated in CASA using the clean task with the högbom algorithm (Högbom 1974) and a natural weighting. Considering the large number of channels in the data cubes (1400), we clean the line cubes using a circular mask spanning the extent of the input sky model to make the imaging process run faster. The images at 800 and 2000 pc have synthesised beam sizes of 0.09 × 0.06 , PA=88 • and 0.07 × 0.06 , PA=90 • , respectively. The reason for the slight increase in the size of the beam for the images at 800 pc is due to the more prominent contribution of side-lobe noise at 800 pc as the source is observed to be more extended. The synthesised beam sizes and root-mean-squared (rms) noise values for both continuum and line observations are summarised in Table 1.

NOEMA
For the synthetic NOEMA observations, we use the uv_fmodel task in the GILDAS 3 software package developed by IRAM and Observatoire de Grenoble. To do so, we create a uv table from the input sky model by providing a uv table that specifies the desired sampling of the antennas in the uv-plane. To be able to make direct comparisons to our NOEMA large program, CORE, we use the uv-coverage of our NOEMA observations in the A, B, and D arrays for source W3(H 2 O) . Six 15-m antennas were used at the time of the observations of W3(H 2 O), with baselines extending the range of 19 to 760 m. The upper right panel of Fig. 2 shows the uv-coverage of the synthetic NOEMA observations. The integration time of the synthetic observations are made in 20-minute intervals for a total time of ∼11 hours spread among the three different arrays. The details of the observing sequence is described in detail in Ahmadi et al. (2018). To achieve the highest possible angular resolution, we cleaned the synthetic visibilities using the MAPPING program of the GILDAS package with the clark algorithm (Clark 1980) and a uniform weighting (robust parameter of 0.1). The images have synthesised beam sizes of 0.44 × 0.34 , PA=44 • . The synthesised beam sizes and rms noise values for both continuum and line observations are summarised in Table 1. The contours for the model simulations are drawn at 1, 5, 9, 20, 40, and 60% of the peak of the emission. The two outer-most contours for the synthetic observations correspond to the 6 and 12σ levels, the rest start at 30σ and increase in steps of 80σ (see Table 1 for the noise values). Synthesised beams are shown in the bottom left corners and scale bars are shown in the bottom right corners of each set of synthetic observations. Note that the intensity scale is different in each panel.

Analysis & Results
In the following, we will present the continuum maps, and study the kinematics, physical conditions, and stability of the disks. We refer to the numerical simulations as 'model simulations' and the synthetic ALMA and NOEMA observations accordingly. Figure 3 shows the 1.3 mm dust continuum models of the simulations at 800 and 2000 pc distance at four different inclinations (10 • , 30 • , 60 • , and 80 • ) and the corresponding synthetic ALMA and NOEMA observations. The synthesised beam sizes, linear resolutions, and rms noise values are summarised in Table 1. At this snapshot, after about 12 kyr of evolution, the fragmented accretion disk spans an extent of ∼1000 au in dust continuum with four fragments seen on scales of ≤500 au, which are actively accreting material from their own smaller accretion disks (described in detail in the following section). Over-dense regions connecting the fragments are a result of gravitational instabilities in the disk which had initially formed spiral arms and resulted in the creation of the fragments. The disk at this stage is still highly dynamic, with the fragments rotating about the central protostar, migrating inwards, and encountering orbital interactions and mass transfer with each other. The center, where the protostar is actively accreting material from the disk and envelope, has a half doughnut shape as a result of the boundary conditions of the numerical simulations. The central sink has a radius of 30 au, the region within which is not included in the computation domain. As a result, the strongest continuum peak is the fragment found to the north-west of the protostar. The arc of material seen in the lower part of the disk shows the ejection of some material from the system. The halo seen around the disk corresponds to the envelope which is much less dense than the disk.

Continuum maps
The synthetic ALMA observations at 800 pc, corresponding to a linear resolution of 60 au, resolve all fragments at all inclinations. The rms noise of the continuum ALMA images at 800 au is higher than the case of the ALMA images at 2000 au as the extended nature of the structure results in more prominent side-lobe noise. Therefore, although the arc of material being ejected in the bottom is visible, it is not detected with 6σ or higher certainty. On the contrary, the synthetic ALMA observations at 2000 pc have lower side-lobe noise and therefore detect this arc with 6σ confidence. At this distance, with a linear resolution of 130 au, ALMA is able to detect all fragments but barely the closest one to the protostar and detangling the fragments at an inclination of 80 • becomes difficult. This disparity in the nature of the noise and the more prominent contribution of side-lobe noise at 800 pc is the reason why the synthesised beam at 800 pc has a slightly larger size than at 2000 pc.
The synthetic NOEMA observations at 800 pc, corresponding to a linear resolution of ∼300 au, are able to resolve all the fragments at 10 • and 30 • inclinations, with only three being resolved at 60 • inclination, and two at a near edge-on view of 80 • inclination. At 2000 pc, corresponding to a linear resolution of ∼800 au, only one structure is resolved regardless of the inclination. The structure however has an elongated shape in the direction of the brightest fragment in the north-west, hinting at the possibility of unresolved fragmentation. At both distances, the disk is seen to have a larger extent due to beam smearing. The angular resolution of NOEMA is expected to improve by roughly a factor of two with the expected baseline extension in the coming years. With the upgraded NOEMA, our synthetic NOEMA observations at 2000 pc will roughly resemble the presented observations at 800 pc.

Kinematics
In the following, we study the kinematics of the disk through moment analysis and position-velocity diagrams.

Moment maps
We create moment maps of CH 3 CN (12 4 − 11 4 ) to study the disk kinematics in more detail. For the model simulations, we add arbitrary Gaussian noise to the cubes and create the moment maps with constraints on the minimum emission level set at 6σ. While there is no need to add noise to the model simulations, it makes for cleaner presentation of the data and easier comparison with maps of the synthetic observations, as well as allowing for the automation of the procedures with the observations. For the synthetic observations, due to the non-Gaussian nature of noise and the large contribution of side-lobe noise in the imaged cubes, we take a higher threshold of 18σ to constrain our study of the kinematics to the disk region. The range of velocities over which the moment maps have been integrated is the same for the model simulations and the synthetic observations. Figure 4 shows the integrated intensity (zeroth moment) map of of CH 3 CN (12 4 − 11 4 ) for the model and synthetic observations, with continuum contour overlays. The distribution of CH 3 CN gas is found to be extended by a few hundred astronomical units beyond the compact dust emission, with the peaks of gas emission coinciding with dust emission peaks as expected. This extension of gas is due to the contribution of the envelope. The envelope itself is rotating and infalling onto a ring-like structure at the location of the centrifugal barrier which can best be seen in the zeroth moment maps of the 10 • and 30 • inclined cubes (further discussed in Sect. 5.2.2). The arc of ejected material in the bottom is also seen in the integrated intensity maps and is well recovered in some of the synthetic observations. The distribution of gas in the synthetic NOEMA observations at 2000 pc are more symmetrically distributed than the dust emission which has an elongated shape. This is because the fraction of gas that is being accreted onto each of the fragments is completely washed out by the large beam and therefore the main contribution left is the larger-scale disk component that is feeding the protostar at the center. Similarly, in the synthetic NOEMA observations at 800 pc the peak of the integrated intensity map is at the center whereas the dust continuum peaks at the location of the fragment to the north-west, since the inner-most 30 au around the accreting protostar is not included in the computation domain, due to the existence of the sink particle, and is therefore not bright in the continuum observations. Figure 5 shows the intensity-weighted peak velocity (first moment) map of CH 3 CN (12 4 − 11 4 ) for the numerical simulations and synthetic observations, with continuum contour overlays. The velocity gradient seen at 10 • inclinations mainly show the contribution from the infalling and rotating envelope with the edges of the envelope shell clearly visible as it has the highest velocity offset from the local standard of rest (LSR) velocity. As the structure gets more inclined, the disk contribution becomes more clearly visible with larger velocity variation seen between the redshifted and blueshifted gas. The synthetic observations are able to recover the kinematics of the gas well, with the envelope and disk contributions becoming heavily blended in the poorly-resolved NOEMA observations making the amplitude of   the velocity gradient smaller. The Yin-Yang appearance of the first moment map is due to the fragments accreting some of the infalling gas which would have otherwise accreted onto the central protostar had the fragments not existed. This can be better seen in Fig. 6 which shows zoom panels of the first moment map of CH 3 CN (12 4 − 11 4 ) on the positions of the fragments for the model simulations at 800 pc and inclined 60 • . The lineof-sight velocity gradient across the central protostar is in the east-west direction and has the largest amplitude as the accretion rate is highest onto the central protostar. The rotation axis of small accretion disks surrounding each of the fragments is not necessarily inherited from the east-west rotational motion of the large-scale disk, but is set according to local dynamics. In the more inclined views towards edge-on, the Yin-Yang effect is less prominent than in the more face-on views as the fragments and their small disks are seen in one plane and start to obscure each other.

Position-velocity diagrams
To study the kinematics in more detail, we create positionvelocity (PV) plots for a cut in the east-west direction, perpendicular to the direction of the rotation axis. Figure 7 shows the PV plots of CH 3 CN (12 4 − 11 4 ) for the model simulations and synthetic observation. The PV plots of the model simulations show the effect of inclination clearly in that the inner disk contribution (r < 250 au) is not strong in the more face-on views. As the source is more inclined towards edge-on, the disk component becomes prominent, with high-velocity contributions seen very close to the central protostar. The inner disk is in clear differential rotation with a Keplerian profile v(r) = √ GM/r about a central object with a mass of ∼10 M as depicted by the green curve. The contribution seen at r > 250 au corresponds to the region where the rotating envelope is falling onto the inner disk which has a lower radial velocity than the infall velocity of the envelope. In recent years, this transition has been observed around low-mass protostars (e.g., Sakai et al. 2014;Oya et al. 2016;Lee et al. 2017;Alves et al. 2017). These studies have found that  Synthesised beams are shown in the bottom left corners of each set of synthetic observations. The contours correspond to 1.37 mm continuum as shown in Fig. 3. The scale of the plots has been chosen for ease of comparison to the 'true' disk mid-plane temperature as shown in Fig. 9. at this position, accretion shocks lead to local heating, causing a change in the chemistry as material gets liberated off the ice mantles to the gas phase. Furthermore, evidence for the centrifugal barrier around a high-mass protostellar envelope has recently been reported by Csengeri et al. (2018). This transition region is estimated to be an order of magnitude closer to the central protostar in the low-mass regime (e.g., 30-50 au for IRAS 16293-2422 Source B: Oya et al. 2018) than in the reported high-mass case (300-800 au: Csengeri et al. 2018). Our estimate of the centrifugal barrier located at r ∼250 au is in agreement with the lower end of the estimates from the observations of Csengeri et al. (2018).
The sharpness of this discontinuity feature seen in the PV plots is a function of resolution and inclination angle. It can be seen in the NOEMA observations at 2000 au that the 'kink' seen in the other PV plots as a result of this transition is completely smoothed out. Furthermore, the velocity profile of the infalling rotating envelope would look different depending on the offset between the line-of-sight and the source position as depicted in Figure 3b of Sakai et al. (2014), which can have an effect on the sharpness of this discontinuity. The PV plots presented in this work have been made for a cut that goes through the protostar as we have knowledge of the precise location of the protostar.
The PV plots of the synthetic ALMA observations resemble the PV plots of the model simulations but smoothed in the position direction as expected. The high-velocity components are still clearly detected in the more inclined cases, but the spreading of emission would mean that the Keplerian curve corresponding to rotation of matter about a ∼10 M star does not perfectly fit anymore, especially in the 2000 pc case. As we move to the synthetic NOEMA observations and the emission gets even more smeared in the position direction, the high-velocity components seen in the more inclined views become hard to detect and the Keplerian curve does not fit well to the observations. This is best seen in the synthetic NOEMA observations at 2000 pc where the inner disk and envelope contributions are completely blended and the resulting PV plots better resemble rigid-body rotation than differential rotation, similar to many observational data in the last decade at comparable resolution.
As temperature is expected to be higher in regions of the disk closer to the protostar, one would expect line transitions which require higher excitation energies to be present in regions closer to the star. One signpost of differential rotation in observations is often seen in the steepening of the slope of the emission in the PV diagrams for the higher excited lines, such that a transition that is tracing the innermost regions of the disk would have higher velocities at smaller radii than the transitions that trace the entire disk region. This has been observed in the case of higher-K transitions of CH 3 CN (Johnston et al. 2015), or vibrationally excited CH 3 CN transitions (Cesaroni et al. 2017). Comparing the PV plots of CH 3 CN (12 2 − 11 2 ) to CH 3 CN (12 4 − 11 4 ) and CH 3 CN (12 6 − 11 6 ), we do not see this steepening of the slope with higher-K levels because all lines trace regions with similar optical depths and therefore probe similar regions.

Temperature distribution
To gain a better understanding of the disk structure, we calculate the rotational temperature and column density of the dense-gas tracer CH 3 CN under the assumption of Local Thermodynamical Equilibrium (LTE). To do so, we use the eXtended CASA Line Analysis Software Suite (XCLASS 4 , Möller et al. 2017) software in CASA. In Fig. 8 we show the resulting rotational temperature maps obtained from pixel-by-pixel XCLASS modelling of CH 3 CN (12 K − 11 K ) K = 4 − 6 for the synthetic observations. The recipe is as outlined in Appendix B of Ahmadi et al. (2018), where we discuss the reasoning behind not including the lower K-ladder transitions in the fitting routine. In these runs, we have also allowed for the size of the source to be a fitting parameter rather than fixing it to a value larger than the beam. This strategy has been tested and yields smoother temperature maps. For each synthetic observation, we calculate the rms noise in the emission-free region in the channel that has the peak of emission for CH 3 CN (12 4 −11 4 ), and only fit pixels with signal higher than 6 times this noise. We note that the spectra at the positions where the envelope is directly falling onto the disk are triply peaked corresponding to red-and blue-shifted infall from the envelope onto the disk with a third component corresponding to the rotational motion of the disk itself. This is particularly prominent in the views closer to face-on. Since performing a three-component fit to the entire map was not feasible, we smoothed the velocity resolution by a factor of 8, reaching a spectral resolution of 2.8 km s −1 , to blend all three components of the lines enough to perform a one-component fit to the spectrally smoothed data cubes. Figure 8 shows the resulting rotational temperature maps for the synthetic observations. Under the LTE assumption, the rotational temperature of CH 3 CN can be accepted as the kinetic temperature of the gas. The median gas temperatures obtained from the synthetic observations at each inclination are listed in Table 2. On average, temperature decreases as the system becomes more inclined, which is expected as the observed column densities are higher in the more inclined views. Also, the arc of material being ejected from the system seen to the bottom of the disk has a higher temperature than its surroundings. For comparison, the face-on view of the disk mid-plane temperature from the model simulations is presented in Fig. 9 on the same logarithmic scale. Mid-plane temperatures at the location of fragments 4 https://xclass.astro.uni-koeln.de can reach values higher than 1000 K, but on average the median temperature for the inner disk (r < 400 au) is 177 K, about 50 K warmer than the outer envelope (400 au < r < 800 au), which has a median temperature of 128 K. In Table 2, we also list median temperature values calculated for the regions within 12σ dust continuum contours (the second outermost contours in Fig. 8) which better correspond to the disk. On average, the temperatures are as expected higher in the disk than in the outer envelope. Moreover, the temperature decreases as angular resolution worsens when the disk structure is no longer recovered. Interestingly, the temperature distributions for the synthetic NOEMA observations at 2000 pc at 10 • and 30 • inclinations are elongated in the same direction as the continuum, while this elongation feature is not seen in the respective zeroth moment maps of CH 3 CN (12 4 − 11 4 ) (bottom right panel of Fig. 4).

Mass estimates
At the snapshot of the numerical simulations under investigation, the protostar has gained a mass of ∼10 M . In the following, we use different techniques often used in observations for determining core and protostellar masses.

Masses from dust emission
Because dust opacity decreases with increasing wavelength, thermal emission of dust at longer wavelengths in the Rayleigh-Jeans limit can trace large column densities and provide estimates for the dust mass if the dust emissivity is known (Hildebrand 1983). Assuming optically thin dust emission at 1.3 mm, for a gas-to-dust mass ratio R, and a dust absorption coefficient κ ν , we can convert the flux density F ν of the continuum observations to a mass via Article number, page 13 of 20 A&A proofs: manuscript no. sims_arXiv_v2 Notes. (a) Median temperature, calculated from the maps presented in Fig. 8. (b) Median temperature within 12σ dust continuum (second outermost) contours as presented in Fig. 8. (c) Mass is calculated from the dust emission using Eq. 5.
where d is the distance to the source, and B ν (T D ) is the Planck function which at the wavelengths under investigation follows the Rayleigh-Jeans law and is linearly dependent on the dust temperature T D . We adopt a value of 150 for the gas-to-dust mass ratio (Draine 2011) and κ ν = 0.9 cm 2 g −1 corresponding to thin ice mantles after 10 5 years of coagulation at a density of 10 6 cm −3 (Ossenkopf & Henning 1994). We assume the kinetic temperature of the gas, derived from the radiative transfer modelling of CH 3 CN (12 K − 11 K ) K = 4 − 6 as presented in Sect. 5.3, to be coupled to the dust temperature. Therefore, we use the temperature maps presented in Fig. 8 together with the continuum maps shown in Fig. 3 to create mass maps. Summing over all pixels, we then obtain a total mass for each synthetic observation, summarised in Table 2. The values are on average between 2.5-3.5 M . For comparison, the total mass within the inner 8400 by 8400 au box for the model simulations is ∼14 M calculated in the same manner by using the face-on continuum map and the mid-plane temperature shown in Fig. 9, which accounts for the total gas mass in this area. The bulk of the remaining mass is on larger scales not included in our analysis. Therefore, the mass estimates from the synthetic observations recover roughly 20% of the total mass as calculated from the model simulations and the bulk of the flux is filtered out by the interferometer.
Considering only the synthetic ALMA observations, the mass estimates at 2000 pc are on average smaller than at 800 pc, by roughly 10%. Since the temperature distribution on average does not vary significantly between the two cases, one would expect the mass estimate to also stay consistent. This small excess of mass corresponds to the contribution from the halo of material found outside of the 6σ continuum contours in the middle panels of Fig. 3. Since the continuum maps at 800 pc are noisier than at 2000 pc (see Table 1) the mass estimates at 800 pc are higher than at 2000 pc. Conversely, the mass estimates at 2000 pc for the synthetic NOEMA observations are roughly 20% higher than at 800 pc. While the average temperature is not significantly different in the two cases (see Table 2), the distribution of temperature in the inner-most regions where the bulk of continuum emission is concentrated is significantly cooler in the 2000 pc case than at 800 pc (see bottom row of Fig. 8), resulting in a higher mass estimate at the farther distance.

Protostellar masses from PV plots
Protostellar masses are often estimated from fitting Keplerian profiles with v(r) = √ GM/r to PV diagrams. In Fig. 7 we show the PV plots for the model simulations and synthetic observations for a cut along the east-west direction going through the center of the sink particle where the protostar is located. While the green curves correspond to the region within which emission is expected if the gas is in a disk in Keplerian rotation about a 10 M star, they do not represents fits to the PV plots. In order to estimate protostellar masses from these PV plots, we use the KeplerFit package 5 (described in Bosco et al. 2019) which uses the method presented in Seifried et al. (2016). The method first identifies the two quadrants of the PV plot that contain most of the emission, and then starting from the positions closest to the protostar moving outwards it determines the highest velocity at which the signal becomes higher than a given threshold. This then determines the boundary which these authors call the 'upper edge' of the PV plot which is then fit with a Keplerian profile to determine the enclosed mass. This method has proven to yield more accurate mass estimates than fitting the pixels with the maximum intensity for each radius or for each velocity channel (Seifried et al. 2016).
It can clearly be seen from the PV plots shown in Fig. 7 that the rotating and infalling envelope (regions outside of the vertical dashed lines) shows different kinematic signatures than the inner disk as discussed in Sect. 5.2.2. More specifically, there exists a discontinuity in the resolved PV plots as a result of the infalling and rotating envelope having a higher infall velocity than the inner disk. Therefore fits to the inner regions would be expected to yield lower mass estimates than to the outer envelope which at a given position has higher velocities. For this reason, we fit the inner and outer part of the disk separately and compare the results to the case where the entire PV diagram is fitted. An example is shown in Fig. 10 for the model simulations at 800 pc. We always fit the 6σ outer edge for a range of radii and always exclude in the fitting the regions closest to the protostar where the velocity decreases with decreasing distance to the protostar as a result of finite resolution (see Appendix A.3.2 of Bosco et al. 2019). The resulting mass estimates from fitting the inner, outer, and entire PV curves for the model simulations and synthetic observations are presented in Fig. 11. Knowing that the central protostar has gained a mass of ∼10 M and that the disk mass is ∼8 M at this snapshot of the simulation, we can test the accuracy of estimating the enclosed masses using this method in each scenario. We will first focus on the general findings.
As rotational velocity is expected to scale with inclination angle i by a factor of sin(i), the mass estimates would be expected to scale by a factor of sin 2 (i). Therefore, it is expected that fitting the kinematics of the more inclined disks would result in higher mass estimates than their less inclined counterparts (grey curves in Fig. 11). While this general trend is seen in our fitting results, we find that the fitted masses are on average higher than the theoretical expectation, similar to the findings of Seifried et al. (2016). For example, one would expect the enclosed mass estimate at 30 • to be 0.25 times lower than the true value of 18 M at 90 • inclination; however, the mass estimate from the fitting the inner region of the model simulation is larger than the expected 4.5 M by 35%. Seifried et al. (2016) attribute this increase to the existence of considerable radial motions which in the more inclined views noticeably contribute to the line-ofsight velocities. It is worth noting that while the 10 M protostar is located at the center of the computational domain at zero offset in the PV diagrams, the disk which roughly contains 8 M spans a range of radii, such that the mass-velocity relation is v(r) = √ G(M * + M disk (r)/r. Therefore, the enclosed mass is not exactly 18 M at all radii and the true theoretical curve lies somewhere between the solid and dashed lines in Fig. 11. Interestingly, the fits to the inner disk of the model simulations inclined to 60 • and 80 • provide mass estimates in this range.
For the model simulations where we fit the PV curves of the model simulations and the ALMA observations at 800 pc, the protostellar mass estimates obtained from fitting the inner disk regions are as expected smaller than the estimates from fitting the outer or entire regions, and provide estimates closest to the expected values. For these cases, fits to the entire region yield masses that are less than the mass estimates obtained from fitting the outer regions. In these cases, fits to the outer regions overestimate the masses because the component that is fitted is mainly the envelope within which more mass is contained and which does not necessarily have a Keplerian-like rotation profile.
For all observations, the method overestimates the mass, even when fitting the inner disk regions. The mass is much more overestimated for disks with lower inclinations (i.e., views closer to face-on). The fact that we generally overestimate the masses especially at lower inclinations can be seen as favourable since without the knowledge of the inclination in most observations, one typically adopts the mass estimate obtained from fitting the PV diagram as the protostellar mass without any corrections for inclination. Therefore, it is somewhat reassuring that fits to the entire region for the model simulations and all ALMA observations, which properly resolve the disk, are quite close to the expected enclosed mass of ∼18 M at all inclinations.
As resolution worsens, in the case of the synthetic ALMA observations at 2000 pc and synthetic NOEMA observations at both distances, the fits to the outer regions provide the lowest mass estimates which are also closer to the expected value. For the ALMA observations, the variations in the mass estimates from fits to different regions is not as significant as for the NOEMA observations at 800 and 2000 pc where the fits to the inner regions can yield masses that are larger than the fits to the outer regions by almost a factor of 2 and 4, respectively.
While the Seifried et al. (2016) method of fitting the outer edge of the PV diagram for determining the enclosed mass is the most accurate method currently used, we find that fitting the entire PV diagram of poorly-resolved observations yields masses that are highly overestimated. This is expected as the emission is smeared over a larger area with the envelope and disk components completely blended, resulting in the PV diagrams not having a Keplerian-like shape (see bottom row of Fig. 7). In such cases, we find that fitting the outer regions of the PV diagram (in our cases, regions beyond 1600 au) provides better estimates for the protostellar mass. The same analysis was done for the PV plots of CH 3 CN (12 K − 11 K ), K = 3, 5, 6 lines and the variations in the mass estimates were marginal as compared with those presented in Fig 11 for the K = 4 transition. The only noticeable difference was that the mass estimates obtained from fitting the PV diagrams of the K = 3 transition were slightly higher than those obtained from the other lines, as this transition is more easily excited than the others and the emission is hence distributed over a larger area. The triangles pointing to the left and right correspond to fits to the inner and outer regions, respectively. The diamonds correspond to fits to the entire PV curves. The transparency of the markers correspond to the χ 2 value of the fit, weighted by the χ 2 value of the best-fit curve for a given region. The more opaque the marker, the better the fit. The solid and dashed grey curves correspond to the true estimate for the mass of the protostar (10 M ), and the mass of the protostar + disk (18 M ), respectively, corrected by sin 2 (i) for each inclination i.

Toomre stability
As mentioned in the introduction, the Toomre Q parameter can be used to study the stability of a disk against gravitational collapse (Toomre 1964). The stabilising force of gas pressure is included in the expression of Q (Eq. 1) via the sound speed, where γ is the adiabatic index with a value of 7/5 for diatomic gas, k B is the Boltzmann constant, µ is the mean molecular weight with a value of 2.8, and m H is the mass of the hydrogen atom. Furthermore, we account for the stabilising force of shear by assuming the disk is in Keplerian-like rotation such that the angular velocity of the disk at a given radius r would be In Fig. 12 we show the 'true' Q map obtained from the model simulations. As expected, we see low Q values at the positions of the fragments since these regions have already collapsed to form protostars. There also exist low Q values in a ring coincident with the edge of the disk. Since the mid-plane temperature is smoothly decreasing towards the edges of the disk (see Fig. 9), one would expect a decrease in the Q parameter in the outskirts of the disk, but the existence of low Q values in a ring-like structure is due to enhanced disk surface density in a ring of material in this region (see Fig. 1). It is important to note that the Toomre analysis is meaningless outside of the disk region, therefore the transition from stable (high Q) on envelope scales to unstable (low Q) at larger radii in the outskirts of the maps has no physical meaning. This transition is best seen in the view closest to edge-on in Fig. 12 at a declination offset of about ±800 au.
In the following, we estimate the distribution of Q for the synthetic observations and start the analysis by adopting the known protostellar mass at this snapshot and correcting the angular velocities for the known inclinations. The effects of not correcting the velocities properly and uncertainties in mass estimates are further discussed in Sections 5.5.1 and 5.5.2, respectively. We assume that the LTE assumption is valid in such a high-density environment, and use the rotational temperature of CH 3 CN, modelled in Sect. 5.3, as representative of the kinetic temperature of the gas. To account for the self-gravity of the disk, we include the mass of the disk within a given radius, M disk (r), to the mass of the protostar M * assumed to be 10 M . The angular velocity is further corrected to account for the inclination of the In particular, the radius at each pixel in a face-on disk is calculated as r = ∆x 2 + ∆y 2 where ∆x is the distance from the protostar in right ascension and ∆y is the distance from the protostar in declination. Since we do not rotate the structure but only incline it about the x-axis, the projected radius is r proj = ∆x 2 + ∆y 2 cos 2 (i).
The mass of the disk is calculated as discussed in Sect. 5.4.1 from the continuum emission maps, by summing up the mass contained within each projected radius using the aperture photometry tools of the photutils package (Bradley et al. 2016) within the astropy package (Astropy Collaboration et al. 2013) in Python. Furthermore, the beam-averaged surface density of the disk is calculated via where S ν is the peak intensity and Ω B is the beam solid angle. In Fig. 13, we show the Q maps for the synthetic observations. Again, as expected, we find low Q values below the critical value of 1 at the positions of the fragments in the synthetic observations which resolve them. This is attributed to the enhanced surface density at the positions of the fragments which have already collapsed to form stars. We also find low Q values (∼1-2) in the arms connecting the fragments as well as the arc of material to the south being expelled from the system. The arms connecting the fragments can be locations where further disk fragmentation can take place. The extremely low Q values above and below where the continuum contours are drawn in the view closest to edge-on (80 • ) are a result of low angular velocities at these positions, since the disk is only extending a small distance in declination. The gas contribution at these locations come from the rotating and infalling envelope and the Toomre analysis is actually not applicable for this envelope medium.
The ring of low Q values seen in the model simulations (see Fig. 12) is not seen in the Toomre maps of the synthetic observations. Since the produced dust continuum images are essentially a convolution of the temperature with the density, the radial decrease of temperature in the outskirts (see Fig. 9) where the ring exists in the density (see Fig. 1) results in this component not being prominent in the final continuum images (see Fig. 3). Therefore, because we use the continuum maps to calculate the surface density of the disks in the Toomre equation, this component is also non-existent. Interestingly, we find high Q values in a ring just outside of the lowest continuum contours, best seen in the synthetic ALMA observations inclined 10 • and 30 • . This is as a result of high gas temperatures in these regions, due to shock heating of gas falling from the rotating envelope onto the disk.
The synthetic NOEMA observations at 800 pc still show low Q values at the positions of the fragments, but the low values expected to be seen in the arms connecting the fragments are absent because the temperature distribution is uniform in the inner regions (see bottom-left panel of Fig. 8). Most interestingly, although the synthetic NOEMA observations at 2000 pc do not resolve the individual fragments, we find low Q values ∼2 in the outskirts of the disk. This is as a result of the structure having a non-circular elongated shape in the direction of the brightest fragment in the north-west in both the continuum and temperature distributions.
The critical value of Q is often taken to be 1 in theoretical calculations, and lowers to ∼0.7 for an isothermal disk of finite thickness (Goldreich & Lynden-Bell 1965;Gammie 2001). Considering the fact that our modelling of the level populations of the CH 3 CN K-ladder provides gas temperatures that may probe layers above the disk mid-plane, we assume that regions where Q is less than ∼2 are unstable against gravitational collapse. The median Q values are listed in Table 2 and follow an opposite trend as the mass estimates. This is because the Toomre parameter is inversely proportional to the surface density of the disk which is calculated in a similar manner as the mass estimates but without the distance dependence.
At first glance, one sees that on average, the Q parameter gets lower with increasing inclination. This is naively expected as in the edge-on views one looks through a larger column. To investigate the effect of inclination further in the resolved observations, the histogram of Q values for synthetic ALMA observations at 2000 pc is shown as an example in the top panels of Fig. 14. We attempt to remove as much contribution as possible from the envelope by only plotting the histogram of Q values for the pixels that lie within 6σ continuum contours which are the outer-most contours drawn on the maps in Fig. 13. The Q distributions in these inner regions do not vary much for the 10 • , 30 • , and 60 • views, but still a large portion of the contribution for the 80 • inclined map comes from the envelope since the disk contribution is just a narrow region spanning a few pixels in declination. The angular velocities in the envelope for the case closest to edge-on are low and therefore result in low Q values there. For comparison, in the bottom panel of Fig. 14 we show the histogram of Q values for the inclined models where only pixels within a radius of 1600 au are included in order to exclude the meaningless contribution from the larger-scale cloud. While the distribution is broader as there are more pixels with low Q values for each fragment, the findings remain the same.  Fig. 13), including only pixels that lie within 6σ continuum contours to show the Q distribution mostly associated with the disk rather than the envelope. Bottom: Histogram of Q values for model simulations (see Fig. 12), including only pixels within a radius of 1600 au.

Inclination of the angular velocity field
In true observations, one seldom has information about the inclination of the source. Sometimes the outflow geometry can help, but high-mass stars form in highly clustered environments, and with the outflow angles widening with time (Beuther & Shepherd 2005), there is a high chance that outflows emanating from different sources may overlap, making the deduction of an inclination angle for the disk difficult. To check how the lack of knowledge about the inclination of the disk would affect our analysis, we have produced an angular velocity map where at each pixel r = ∆x 2 + ∆y 2 , which essentially means we assume the angular velocity at each position is what it would be if the disk was edge-on in that direction. Figure 15 shows the Q maps for the synthetic observations without the inclination correction. In the resolved ALMA observations at both distances and NOEMA observations at 800 pc, we still find low Q values as expected at the positions of the fragments, but the higher inclined views show higher Q values at the locations of the fragments than in the inclination-corrected cases, precisely because the angular velocities at a given point above and below the zero declination line are higher in the uncorrected maps. Interestingly, in the synthetic NOEMA observations at 2000 pc where the fragments are not resolved, we still see the low Q values in the outskirts of the structure. Therefore, the conclusion from the corrected case still holds that we are able to predict fragmentation of the disk without actually resolving the fragments.

Uncertainties in the mass of the central object
In our Toomre analysis, we have assumed the mass of the protostar at the center to be 10 M as we have access to this information. In Sect. 5.4.2 we outlined how the mass of the central object is often estimated from the kinematics of the disk, and explained that with loss of resolution, the shape of the rotation curve of the disk changes from differential to a more rigid-bodylike appearance. Therefore, fitting Keplerian-like curves to the PV plots of poorly-resolved structures would overestimate the Fig. 15. Same as Fig. 13 but assuming the inclination of the object is not known, therefore, with the velocities not corrected to the actual rotation velocities. mass of the central object. Since the mass of the protostar has a square root dependence in the expression of Q (Eq. 1), overestimating the protostellar mass by a factor of 4 would yield Q values that are on average higher by a factor of two, making the disks seem more stable.

Summary and Conclusions
In this work, we have presented synthetic ALMA and NOEMA observations for a high resolution 3D radiation-hydrodynamic simulation of a high-mass protostellar disk that fragments into a highly dynamic system. At the snapshot of the simulation at 12 kyr, the system comprises of a central protostar and four companion fragments on scales ≤500 au, all accreting material from the infalling and rotating envelope which feeds the disk that is in Keplerian rotation. Smaller accretion disks are formed around each of the fragments, through which some of the large-scale disk and envelope material is accreted onto the fragments. The aim of the investigation has been to study the effect of inclination and spatial resolution in recovering disk structure and to study the stability of the disk in order to predict disk fragmentation and better understand recent and future high-resolution observations. The following is a summary of our findings: -Synthetic ALMA observation at 800 pc are able to detect all fragments at all inclinations in both continuum and line emission at 1.3 mm. At 2000 pc, the ALMA observations at 80 • inclination cannot resolve the closest fragment to the protostar. The NOEMA observations at 800 pc can resolve all fragments at 10 • and 30 • inclinations but not all at higher inclinations. At 2000 pc, NOEMA observations only show a single structure, which looks slightly elongated towards the brightest fragment. -In the zeroth moment maps of CH 3 CN (12-11) rotational transitions, the peak of emission is at the position of the protostar. For the ALMA observations at both 800 and 2000 pc, we are able to detect the rotating envelope which is infalling onto a ring-like structure at the position of the centrifugal barrier. The effect is best seen in the 10 • and 30 • inclined views. -In the first moment maps of CH 3 CN (12-11) rotational transitions, the disk contribution is best seen in the higher in-clined views as expected. There exist small accretion disks around each fragment feeding them part of the disk material; this leads to the large-scale velocity field of the disk having a Yin-Yang shape. Envelope and disk contributions are heavily blended in the poorly resolved case of NOEMA observations at 2000 pc (linear resolutions ∼800 au) where the emission is also smeared over a larger distance, making the amplitude of the observed velocity gradient smaller. -PV diagrams of CH 3 CN (12-11) rotational transitions for a cut through the disk perpendicular to its rotation axis show the differential (Keplerian-like) rotation of material best in the close to edge-on view. At a radius of ∼250 au a discontinuity is seen, corresponding to the position of the centrifugal barrier because the infalling envelope has a higher infalling velocity than the radial velocity of the disk. As the entire structure becomes less resolved in the NOEMA observations, the high-velocity components of the disk become washed out and blended with the envelope component, making the PV diagrams resemble rigid-body-like rotation. -We fit the emission of CH 3 CN (12 K − 11 K ) K = 4 − 6 with XCLASS to obtain rotational temperature maps for the disks. On average we find temperatures in the range of 100-300 K, decreasing as the system becomes more inclined. The envelope is roughly a factor of 2 cooler than the disk. In the poorly resolved NOEMA observations, the envelope and inner disk contributions become blended and the difference between the temperature of the inner and outer regions is not as apparent as in the resolved cases. -Assuming gas and dust temperatures are coupled, mass estimates from dust for the synthetic observations account for roughly 20% of the mass calculated from the model simulations for the inner 8400 by 8400 au and the rest is filtered out by the interferometers. -We make use of the method presented in Seifried et al. (2016) to fit the PV plots with Keplerian models to obtain protostellar mass estimates. Since there exists a discontinuity in the 6σ outer edge of the PV plots as a result of material falling from the envelope onto the disk, we fit the inner disk and outer envelope regions separately. In the case of wellresolved disks with Keplerian-like profiles, we find the fits to the inner regions the most accurate and the fits to the entire PV curves adequate for estimating enclosed masses. For poorly-resolved disks, the fits to the outer regions of the PV diagrams yield much more accurate mass estimates as opposed to fits to the entire region which overestimate the mass. The factor by which the mass is overestimated depends on how poorly the disk is resolved. -Studying the Toomre stability of the disks, we find low Q values below the critical value for stability against gravitational collapse at the positions of the fragments and in the arms connecting the fragments for the resolved synthetic ALMA observations. For the synthetic NOEMA observations at 800 pc (linear resolutions ∼300 au), low Q values are seen at the locations of the fragments but the unstable region in the disk is no longer seen since the temperature structure is uniformly warm for the inner disk. For the NOEMA observations at 2000 pc (linear resolutions ∼800 au) where only one structure is resolved, we find low Q values in the outskirts of the disk. This is an important result as it implies that despite the lack of ability in resolving any of the fragments, we are able to predict that the disk is unstable and fragmenting. -Assuming that the inclination of the disk is not known, recalculating the Q maps using Keplerian velocities of an edge-on disk to all radial directions for all inclinations provide slightly higher values for the Q parameter but the conclusions from the previous point still hold. This is reassuring as disk inclinations are difficult to establish.
In this work, we showcased the potential and limitations to study disks in high-mass star formation with current (mm) interferometers. We benchmarked a method to study the stability of such disks, showing that even with poorly resolved observations it is possible to predict whether a disk is prone to fragmentation. We aim to apply our methods to the study of a large sample of high-mass disk candidates within the IRAM large program CORE in an upcoming paper.