On the ALMA observability of nascent massive multiple systems formed by gravitational instability

TL;DR

This study uses radiative transfer calculations and synthetic ALMA images to demonstrate that substructures like gaseous clumps and spiral arms in massive young stellar objects' accretion discs are observable with current ALMA capabilities, revealing insights into disc fragmentation.

Contribution

The paper provides the first detailed simulation-based analysis of ALMA's ability to detect disc substructures in forming massive stars, highlighting the impact of gravitational instability on observability.

Findings

ALMA can resolve disc substructures at 1.2 mm within 1-2 kpc.

Substructures are observable regardless of viewing angle.

Disc evolution increases the likelihood of detecting clumps.

Abstract

Massive young stellar object (MYSOs) form during the collapse of high-mass pre-stellar cores, where infalling molecular material is accreted through a centrifugally-balanced accretion disc that is subject to efficient gravitational instabilities. In the resulting fragmented accretion disc of the MYSO, gaseous clumps and low-mass stellar companions can form, which will influence the future evolution of massive protostars in the Hertzsprung-Russell diagram. We perform dust continuum radiative transfer calculations and compute synthetic images of disc structures modelled by the gravito-radiation-hydrodynamics simulation of a forming MYSO, in order to investigate the Atacama Large Millimeter/submillimeter Array (ALMA) observability of circumstellar gaseous clumps and forming multiple systems. Both spiral arms and gaseous clumps located at ~a few 100 au from the protostar can be resolved by interferometric ALMA Cycle 7 C43-8 and C43-10 observations at band 6 (1.2 mm), using a maximal 0.015" beam angular resolution and at least 10-30 min exposure time for sources at distances of 1-2 kpc. Our study shows that substructures are observable regardless of their viewing geometry or can be inferred in the case of an edge-viewed disc. The observation probability of the clumps increases with the gradually increasing efficiency of gravitational instability at work as the disc evolves. As a consequence, large discs around MYSOs close to the zero-age-main-sequence line exhibit more substructures than at the end of the gravitational collapse. Our results motivate further observational campaigns devoted to the close surroundings of the massive protostars S255IR-NIRS3 and NGC 6334I-MM1, whose recent outbursts are a probable signature of disc fragmentation and accretion variability.

Tables

Table 1. Table 1: Accretion disc models selected from  Meyer et al. ( 2018 ) . The table gives the simulation time (in kyr kyr \rm kyr ), protostellar mass (in M ⊙ subscript M direct-product \rm M_{\odot} ), the total protostellar luminosity, the protostellar radius and a brief description of the disc structure for each of the considered disc models.

$\mathrm{Model}$	$t$ ($\mathrm{kyr}$)	$M_{\star }$ ($M_{\odot }$)	$L_{\star }$ ($L_{\odot }$)	$R_{\star }$ ($R_{\odot }$)	$T_{\mathrm{eff}}$ ($K$)	description
1	20.1	$8.6$	$283.0$	$22.3$	$5020.0$	young, stable disc
2	23.0	$10.7$	$650.0$	$30.6$	$5275.0$	stable disc with first gravitational perturbation
3	26.4	$12.1$	$8359.0$	$43.6$	$8374.0$	disc with spiral arms and migrating gaseous clump
4	30.7	$16.6$	$27973.0$	$32.0$	$13210.0$	disc with multiple ring-like spiral arms hosting clumps
5	31.1	$16.7$	$28625.0$	$30.6$	$13586.0$	strongly fragmenting disc with multiple clumps
6	32.1	$20.0$	$226454.0$	$389.0$	$6393.0$	extended strongly fragmented disc with filamentary spiral arms and clumps

Table 2. Table 2: Adopted properties of the bursting MYSOs S255IR-NIRS3  ( Caratti o Garatti et al. , 2017 ) and (b) NGC 6334I-MM1  ( Hunter et al. , 2017 ) .

$\mathrm{Name}$	$M_{\star }$ ($M_{\odot }$)	$L_{\star }$ ($L_{\odot }$)	$R_{\star }$ ($R_{\odot }$)	$T_{\mathrm{eff}}$ ($K$)	$d$ ($\mathrm{kpc}$)	$\varphi$	References
S255IR-NIRS3	$20$	$2.9\times {10}^4$	$10$	${23800}^a$	$1.8$	$80°$	Caratti o Garatti et al. (2017), Meyer et al. (2019)
NGC 6334I-MM1	$20$	$1.5\times {10}^{5a}$	$350$	$6100$	$1.3$	${\varphi }^b$	Forgan et al. (2016), Hunter et al. (2017)

Full text

The observation probability of the clumps increases with the gradually increasing efficiency of gravitational instability at work as the disc evolves. As a consequence, large discs around MYSOs close to the zero-age-main-sequence line exhibit more substructures than at the end of the gravitational collapse. Our results motivate further observational campaigns devoted to the close surroundings of the massive protostars S255IR-NIRS3 and NGC 6334I-MM1, whose recent outbursts are a probable signature of disc fragmentation and accretion variability.

1 Introduction

Multiplicity is an indissociable characteristic of massive star evolution. Observational studies of OB stars indicate that multiplicity is an intrinsic feature, which strongly influences the evolution of massive stars in the Hertzsprung-Russell diagram (Sana et al., 2012). A wide dispersion of the semi-major axis in hierarchical massive stellar systems, which span a range from a small fraction of $\rm au$ for close/spectroscopic companions (Mahy et al., 2013; Kobulnicky et al., 2014; Chini et al., 2012) to up to a few $100\,\rm au$ in the context of massive proto-binaries (Kraus et al., 2017), asks the puzzling question of their formation channel. This must be prior to the launching of the strong stellar winds at the zero-age-main-sequence (ZAMS), which dissipate the disc and therefore limit the companion formation timescale (Vink et al., 2000). On the one hand, investigations of the gravitational capture scenario of stellar systems with N-body simulations failed in explaining this high occurence of the small separations observed in spectroscopic companions to massive stars (Fujii & Portegies Zwart, 2011). On the other hand, the disc fragmentation approach by means of multi-dimensional hydrodynamics simulations coupled to Lagrangian sink particles that spot secondary star formation in self-gravitating accretion discs suffers from the dramatic qualitative algorithm-dependence of the results (Klassen et al., 2016). Therefore, understanding the early formation mechanism of hierarchical multiple systems around massive stars by confronting numerical models to observations is of prime interests in the study of the pre-ZAMS (Zinnecker & Yorke, 2007) but also the pre-supernovae evolution of massive stars (Langer, 2012).

State-of-the-art numerical models of the formation of MYSOs predicted the existence of an accretion disc connected with bipolar cavities filled by ionized stellar wind and outflow material (Seifried et al., 2011; Harries, 2015; Harries et al., 2017). The recent gravito-radiation-hydrodynamical simulations of the formation of massive protostars in Meyer et al. (2017) and Meyer et al. (2018) showed that it is possible to obtain self-consistent solutions for multiplicity in the high-mass regime of star formation without using artificial sink particles, as long as the spatial resolution of the inner disc region is sufficiently high. Those high-resolution simulations of $\simeq 100\,\rm M_{\odot}$ collapsing pre-stellar cores also demonstrate that even irradiated, massive, self-gravitating accretion discs around massive young objects are subject to gravitational fragmentation and generate dense spiral arms in which gaseous clumps form. These clumps can contract further and reach the molecular hydrogen dissociation temperature while migrating inward onto the central forming massive object. The migrating clump can experience a second collapse and become low-mass stars, progenitors of the future companions and/or close (spectroscopic) companions to the massive stars.

The series of studies of Meyer et al. (2017) and Meyer et al. (2018) provides a self-consistent picture unifying disc fragmentation and binary formation in the context of massive young stellar objects. The model also predicts that a fraction of the migrating gaseous clump material can directly fall onto the protostellar surface and produce accretion-driven outbursts as monitored in the MYSOs S255IR-NIRS3 and NGC 6334I-MM1. In this picture, accretion bursts are a signature of the presence of a self-gravitating disc in which efficient gravitational instability occurs. The bursts happen in series of multiple eruptions, whose intensity scales with the quantity of accreted material (Meyer et al., 2019). These episodic events generate modifications in the internal protostellar structure leading to excursions in the Hertzsprung-Russel diagram (Meyer et al., 2019). Note that this burst mode of massive star formation has an scaled-up equivalent in the low-mass regime, as outlined in a series of papers of Vorobyov & Basu (2010); Vorobyov et al. (2013); Dong et al. (2016); Elbakyan et al. (2019), as well as in the primordial mass regime of star formation, see Greif et al. (2012) and Vorobyov et al. (2013). The remaining question is therefore, can such self-gravitating effects be directly observed and what can we learn from the disc density structure about the high-mass star formation process ?

Searching for the direct observational signatures of clumpy structures in massive star-forming regions is a very active field, as modern facilities are now able to probe the highly-opaque pre-stellar cores in which MYSOs grow. The last decade saw the report of observations of variabilities in the accretion flow onto massive protostars (Keto & Wood, 2006; Stecklum et al., 2017) and in the pulsed bipolar outflows which release ionized gas (Cunningham et al., 2009; Cesaroni et al., 2010; Caratti o Garatti et al., 2015; Purser et al., 2016; Reiter et al., 2017; Burns et al., 2017; Burns, 2018; Purser et al., 2018; Samal et al., 2018), together with proofs of disc-jet associations such as G023.01-00.41 (Sanna et al., 2018). Direct observation of Keplerian discs around MYSOs (Johnston et al., 2015; Ilee et al., 2016; Forgan et al., 2016; Ginsburg et al., 2018; Maud et al., 2018) and evidence of filamentary spirals feeding the candidate disc W33A MM1-Main (Maud et al., 2017) and of infalling clumps in the systems G350.69-0.49 (Chen et al., 2017) and G11.92-0.61 MM1 (Ilee et al., 2018) revealed the existence of substructures in discs that are similar to those modelled in Meyer et al. (2018). Furthermore, the investigations for molecular emission from collapsing clouds and from the close surroundings of high-mass protostars confirm the possibility of fragmentation and the presence of inhomogeneities in their discs (Beuther et al., 2019; Ahmadi et al., 2018), without excluding the possibility of magnetic self-regulation of the inner disc fragmentation (Beuther et al., 2018) as predicted by analytic and numerical studies (Hennebelle et al., 2016).

This study aims at testing the observability of the overdensities and various substructures such as spiral arms and gaseous clumps that are predicted by Meyer et al. (2017) and Meyer et al. (2018). We adopt the standard approach of radiative transfer calculations against dust opacities, performed on the basis of the accretion disc density and temperature fields from numerical simulations of Meyer et al. (2018) and further post-processed in order to obtain synthetic interferometric Atacama Large Millimeter/submillimeter Array (alma) band 6 ($1.2\,\rm mm$) images. This corresponds to the typical waveband at which the surroundings of MYSOs are observed, as it offers a good compromise between high spatial resolution and phase calibration problems occurring at longer interferometric baselines. This approach has been adopted to predict the observability of gravitational instabilities in discs around low-mass stars (Vorobyov et al., 2013; Dong et al., 2016; Seifried et al., 2016). Only very few similar studies have been conducted in the context of high-mass stars, either using artificial sink particles as tracers of multiplicity in the disc (Krumholz et al., 2007) or on the basis of simplistic, analytic disc models (Jankovic et al., 2019). Our synthetic images predict what disc morphologies should be observable with alma and constitute prognostications that are directly comparable with future high angular resolution observations campaigns.

Our work is organized as follows. In Section 2, we review the methods used to perform numerical hydrodynamical simulations, radiation transfer calculations and to compute synthetic images of the fragmented circumstellar medium of a massive protostar. In Section 3, we present our results as a series of synthetic images of accretion discs one should expect to observe with alma in the surroundings of forming MYSOs. We consider several viewing angles of the objects with respect to the plane of the sky. Our synthetic images are further discussed in Section 4. Finally, we provide our conclusions in Section 5.

2 Method

In the following paragraphs, we remind the reader about the numerical methods used to carry out our hydrodynamics disc models. The structure and accretion properties of the discs are extracted from the simulation of the collapse of a solid-body-rotating pre-stellar core that produces a central, single massive young stellar object (Meyer et al., 2018). The accretion history onto the MYSO is used to feed a stellar evolution code, in order to self-consistently derive its protostellar surface properties. Then, we present the radiation transfer methods utilised to produce dust continuum millimeter images of the discs, by post-processing the density and temperature fields of the hydrodynamical circumstellar disc model. Finally, we review our method for retrieving synthetic alma images, allowing us to transform the radiation transfer calculations into realistic alma synthetic observations and to discuss the observability of the clumps in the disc as a function of, e.g. the protostellar formation phases.

2.1 Hydrodynamical simulations

This study focuses on the high-resolution gravito-radiation-hydrodynamics simulation Run-1-hr of Meyer et al. (2018). It is a numerical model of the gravitational collapse of a $100\,\rm M_{\odot}$ solid-body-rotating cold pre-stellar core of uniform initial temperature $T_{\rm c}=10\,\rm K$ and of density distribution $\rho(r)\propto r^{-3/2}$, with $r$ the radial coordinate. The initial ratio of kinetic-by-gravitational energy of the pre-stellar core was set to $\beta=4\,\%$. The spherical midplane-symmetric computational domain onto which the run has been performed has an inner radius $r_{\rm in}=10\,\rm au$ and an outer radius $R_{\rm c}=0.1\,\rm pc$, respectively. A grid maps the domain $[r_{\rm in},R_{\rm c}]\times[0,\pi/2]\times[0,2\pi]$ with $N_{\rm r}=256\times\,N_{\rm\phi}=41\times\,N_{\rm\theta}=256$ grid cells, expanding logarithmically along the radial direction $r$, going as a cosine in the polar direction $\phi$ and uniformly spaced along the azimuthal direction $\theta$. The inner hole is a semi-permeable sink cell attached onto the origin of the domain whereas outflow boundary conditions are assigned to the outer radius. Therefore, the material $\dot{M}$ that is lost through the sink cell allows us to calculate the accretion rate onto the protostar, while the protostellar properties such as its stellar radius and its photospheric luminosity are time-dependently estimated using the pre-calculated protostellar evolutionary tracks of Hosokawa & Omukai (2009). Because of time-step restrictions due to the high-resolution of the grid, we computed the collapse of the core together with the initial formation phase of the circumstellar disc of the protostar up to $t_{\rm end}=\,32.1\rm kyr$.

We integrate the system by solving the equations of gravito-radiation-hydrodynamics with the pluto code111http://plutocode.ph.unito.it/ (Mignone et al., 2007, 2012), that has been modified to account for (i) the protostellar feedback and (ii) the self-gravity of the gas. The direct proto-stellar irradiation feedback of the central star as well as the inner disc radiative transport are taken into account using the gray approximation which has been adapted from the publicly-available scheme of Kolb et al. (2013), see method section of Meyer et al. (2018). This two-fold algorithm intelligently ray-traces photon packages from the protostellar atmosphere and then diffuses their energy by flux-limited propagation into the accretion disc. However, this scheme allows us to accurately treat both the inner heating and the outer cooling of the irradiated disc surrounding our MYSO (Vaidya et al., 2011). We use a constant gas opacity of $0.01\,\rm cm^{2}\,g^{−1}$ and the tabulated dust opacities of Laor & Draine (1993) in order to account for the attenuation of the radiation field. Gas and dust temperatures are calculated assuming the equilibrium between the temperature of the silicate dust grains and the total protostellar irradiation and disc radiation fields. Any turbulent viscosity is neglected in the computational domain and we assume that the most efficient mechanism for angular momentum transport are the gravitational torques generated in the self-gravitating accretion disc that surrounds the protostar. The gravity of the central star is taken into account by calculating its total gravitational potential and our simulations include the self-gravity of the circumstellar gas by directly solving the Poisson equation. Finally, we refer the reader who is interested in further details about the numerical method to Meyer et al. (2018).

2.2 Stellar evolution calculations

The modelled accretion rate history is used as initial conditions to compute the evolutionary tracks in the Hertzsprung-Russel diagram of the corresponding episodically-accreting massive protostar, following the method of Meyer et al. (2019). To this end, a one-dimensional stellar evolution calculation is carried out with the hydrostatic genec (i.e. geneva code) code (Eggenberger et al., 2008), which was updated with respect to disc accretion physics for the study of pre-main-sequence high-mass stars (Haemmerlé, 2014; Haemmerlé et al., 2016). The accretion mechanism is therein treated within the so-called cold disc accretion approximation (Palla & Stahler, 1992). The inner disc region is assumed to be geometrically thin as the accreted material falls onto the surface of the protostar and the implementation of Haemmerlé (2014) showed full consistency with the original results of Hosokawa et al. (2010). Any entropy excess is immediately radiated away and the circumstellar material is advected inside the protostar, i.e. the simulation considers that the thermal properties of the inner edge of the disc are similar to those of the suface layer of the MYSOs. This approach is the lower limit on the entropy attained by the protostar throughout the accretion mechanism, the upper limit being the so-called spherical (or hot) accretion scenario (Hosokawa et al., 2010). The calculation of the stellar structure is performed with the Henyey method, using the Lagrangian formulation (Haemmerlé, 2014) at solar metallicity (Z=0.014) and making use of the abundances of Asplund et al. (2005) and Cunha et al. (2006), together with the deuterium mass fractions described in the studies of Norberg & Maeder (2000) and Behrend & Maeder (2001). The stellar structure calculation is started with a fully convective stellar embryo (Haemmerlé & Peters, 2016) and make use of overshooting and of the Schwarzschild criterion in the treatment of the internal stellar convection. The outcomes of the calculation are evolution of $L_{\odot}$, $R_{\star}$ and $T_{\rm eff}$ as a function of time.

2.3 Radiative transfer calculations

A series of characteristic simulation snapshots is selected from the simulation Run-1-hr of Meyer et al. (2018) and they are post-processed by performing radiative transfer calculations against dust opacity with the code radmc-3d222http://www.ita.uni-heidelberg.de/ dullemond/software/radmc-3d/ (Dullemond, 2012). We make use of the Laor & Draine (1993) dust mixture based of silicates crystals. Hence, our simulated emission is calculated with the same dust opacity that was used in the gravito-radiation-hydrodynamical simulations. In our approach, radmc-3d directly uses dust density and temperature fields imported from the pluto outputs and we do not perform any Monte-Carlo simulation prior to the image building, see comparison of both methods in Section 4.1. The interface between the pluto and radmc-3d codes is an interpolation between the above described non-uniform mesh of pluto and the uniform spherically-symmetric grid of radmc-3d which reconstructs the full three-dimensional structure of the disc taking into account the midplane-symmetry of the hydrodynamical models. The protostar is located at the origin of both grids and the radmc-3d calculations is then effectuated concentrating onto the inner $1000$-$3000\,\rm au$ regions of the disc, such that the corresponding map has $2000\times 2000$ cells along the horizontal and vertical directions of the sky plane, respectively.

The accretion discs are irradiated by packages of $5\times 10^{6}$ photons ray-traced from the protostellar atmosphere of radius $r_{\rm in}$ to the outer disc region, at a radius of a few $1000\,\rm au$. The protostar is assumed to be a black body radiator of effective temperature $T_{\rm eff}$, bolometric luminosity $L_{\star}$ and stellar radius $R_{\star}$ with values taken from the outcomes of the stellar evolution calculation at the time of the selected snapshots. The radiative transfer calculations are performed within the anisotropic scattering approximation using the Henyey-Greenstein formula. For each of the selected simulation snapshots, we build $1.2\,\rm mm$ emission maps by projecting the emitted flux. This is the typical waveband at which observation campaigns for the study of accretion discs around MYSOs are currently conducted, see e.g. Maud et al. (2018), as this waveband offers a good compromise between high angular resolution and atmosphere phase stability. The distance between the source and the observer is set to $1.0\,\rm kpc$, that is of the order of the distance to the closest massive star-forming regions such as Orion ($0.41\,\rm kpc$), Cep A ($0.70\,\rm kpc$) and VY CMa ($1.14\,\rm kpc$), see Table 4 in Reid et al. (2009), as the angular resolution reached by alma when observing them is the highest possible. Additionally, we also investigate the observability of disc substructures for MYSOs located at a distance of $2.0\,\rm kpc$ which corresponds to further afield star-forming regions such as W3(OH), see Section 3.5. For each of the 6 selected simulation shapshots, we use the possibility of radmc-3d to compute $1.2\,\rm mm$ synthetic images for three different system inclination angles, namely face-on ($\phi=0^{\circ}$), intermediate inclination ($\phi=45^{\circ}$), and edge-on ($\phi=90^{\circ}$). The disc models and their emission properties are presented in detail in Section 3.

2.4 Synthetic observables

We further post-process the radmc-3d radiation transfer calculations with the Common Astronomy Software Applications casa333https://casa.nrao.edu/ (McMullin et al., 2007) in order to obtain synthetic fluxes of the disc models as they would look like if seen by alma observations.

The successive technical updates and enhancements of alma are labelled as ”Cycles” and to each of them corresponds a call for proposals followed by a series of observations. For each cycle, alma operates at different wavebands (the ”bands”). The one of interest in our study is band 6, corresponding to $1.2\,\rm mm$ (Brown et al., 2004; Wootten & Thompson, 2009). Observations are performed with a given spatial configuration of all the antennas constituting the alma interferometer. The more compact the configuration the larger the beam size probing the sky, while the more extended the antennas configuration the better the spatial resolution in the observed images. The antennas are available in two categories, distinguishable by the diameter of their parabola, i.e. $7\,\rm m$ and $12\,\rm m$, respectively. The maximum number of $12\,\rm m$-antenna that can be used simultaneously is 43. We are interested in the $12\,\rm m$ antenna arrays because these 43 antennas can be spread over a wider surface on the Llano de Chajnantor plateau, offering a larger maximal separation between two antennas (the maximum baseline) and hence a smaller beam size to probe the internal structure of our accretion discs. The alma Cycle 7 offers ten $12\,\rm m$ antenna array configurations, labelled from 1 (most compact) to 10 (most extended). We perform synthetic images at two of the most extended configurations (8 and 10) of the 43 different $12\,\rm m$-antenna, also called configurations C43-8 and C43-10. We therefore model images of our discs as alma Cycle 7 configuration 8 and 10 observations at band 6 ($1.2\,\rm mm$), respectively.

The radmc-3d to casa interface is effectuated via the radmc3dPy python package available with the current distribution of radmc-3d. It allows the user to directly output the radiative transfer calculations as standard FITS files, which header can be read by the casa software. We make use the simalma task of casa with a central frequency of the combined continuum emission of $249.827\,\rm Ghz$ ($1.2\,\rm mm$, alma band 6) and a channelwidth of $50.0\,\rm Mhz$.

The different used long-baseline configurations result in angular resolutions of $0.029^{\prime\prime}$ (antenna configuration 8 with maximal baseline of $8.5\,\rm km$) and $0.015^{\prime\prime}$ (antenna configuration 10 with maximal baseline of $16.2\,\rm km$), respectively.

The observation integration time is chosen to be $600\,\rm s$ ($10\,\rm min$) of exposition time, as longer durations do not provide better images (see Section 4). The precipitable water vapor (pwv) parameter is chosen to be $0.5\,\rm mm$, which represent the amount of pwv at the Llano de Chajnantor plateau in the Chilean Andes during $25\,\%$ of the time during which the Atacama Large Millimeter/submillimeter Array (alma) is operated.

3 Results

In this section, we begin by presenting the samples of the hydrodynamical simulation which we post-process to obtain synthetic $1.2\,\rm mm$ interferometric images of dust continuum emission. Then, we discuss the time effects due to the disc evolution and the effects of the disc inclination angle with respect to the plane of sky.

3.1 Hydrodynamical model and stellar properties

The simulation model Run1-hr of Meyer et al. (2018) begins with the gravitational collapse of $100\,\rm M_{\odot}$ of pre-stellar rotating molecular material onto the stellar embryo. At the end of the free-fall collapse phase, the infalling material ends on a centrifugally-balanced disc, from which the gas is subsequently transferred to the growing protostar. Fig. 1 plots the accretion rate history onto the central massive protostar (solid blue line, in $\rm M_{\odot}\,\rm yr^{-1}$) together with the evolution of the protostellar mass (dashed red line, in $\rm M_{\odot}$) that is calculated as the integrated disc-to-star mass transfer rate through the sink cell. The vertical thin black line indicates the onset of disc formation, when the free-fall collapse of the envelope material onto the protostar stops and the star begins to gain its mass exclusively via accretion from its surrounding disc. After the initial infall of material, the collapse of the parent pre-stellar core material generates an initial increase of the mass flux through the inner boundary at a time $\approx 2\,\rm kyr$. The accretion rate then reaches the standard value predicted for MYSOs of $10^{-3}\,\rm M_{\odot}\,\rm yr^{-1}$ (Hosokawa & Omukai, 2009) up to the onset of the disc formation happening at $\approx 12\,\rm kyr$. Variabilities in the accretion flow begin right after the disc formation and the accretion rate history exhibits numerous peaks of growing intensity as the MYSO becomes heavier.

This is not caused by different inner boundary conditions, but simply reflects the time-dependent azimuthal anisotropies induced in the accretion flow by the disc evolution.

The efficient gravitational instabilities produce complex substructures in the disc, such as overdense spiral arms in which gaseous clumps of various morphologies form at radii of $\sim 100\,\rm au$ and inward-migrate down to the central protostar. This produces luminosity outbursts via the mechanism revealed in Meyer et al. (2017) for massive stars. These outbursts are responsible for step-like increases in the stellar mass evolution (see thick dotted red line) due to the fast accretion of dense circumstellar material inside the sink cell. A more detailed description of the evolution of circumstellar discs around young massive stars irradiating their self-gravitating discs can be found in our precedent study (Meyer et al., 2018). In Fig. 1, several magenta dots mark the time instances of the chosen simulation snapshots considered in this study. Note that the young star becomes, by definition, a massive object when $M_{\star}=8\,\rm M_{\odot}$. Hence, our selected disc models are all in the high-mass regime.

Fig. 2 plots the stellar surface properties obtained by post-processing the accretion rate history displayed in Fig. 1 with the genec evolution code, with the evolution of the protostellar iternal photospheric luminosity (a), radius (b), and effective temperature (c) as a function of the stellar age of the MYSO. No notable variations happen up to the end of the free-fall collapse and during the early phase of the disc formation at times $\approx 20\,\rm kyr$. A slight monotonical increase of the photospheric luminosity $L_{\star}$ and of the protostellar radius $R_{\star}$ happens as the deuterieum burning keeps the core temperature almost constant. When the evolution of the accretion rate exhibits sudden variations in response to the accretion of gaseous circumstellar clumps formed in the fragmented accretion disc, it generates changes in the internal structure resulting in the bloating of the protostellar radius (Haemmerlé et al., 2016). At this stage, the entropy of the external layers is high and any episodic deposit of mass on it consequently induces an augmentation of $T_{\rm eff}$ and $L_{\star}$, as described in Meyer et al. (2019). During each strong accretion burst, the MYSOs episodically experiences changes in the photospheric properties by becoming bigger and cooler. Fig. 3 shows the evolutionary track of our MYSO in the Hertzsprung-Russell diagram and compares it with the pre-zero-age-main-sequence (pre-ZAMS) tracks calculated with constant accretion rates (thin solid black lines) of Haemmerlé et al. (2016) and with the ZAMS track (thick dashed black line) of Ekström et al. (2012). The grey solid lines are isoradii. As detailed in Meyer et al. (2019), each strong accretion event reaching rates $\geq 10^{-2}\,\rm M_{\odot}\,\rm yr^{-1}$ induces important rapid changes in the pre-main-sequence evolutionary track of the MYSOs. Following an internal redistribution of entropy in the upper layer of the stellar structure, $R_{\star}$ increases. Successive evolutionary loops toward the red part of the Hertzsprung-Russell diagram occur during these bloating periods. The protostellar properties, required for the radiative transfer calculations, are reported for the 6 selected simulation snapshots in our Table 1.

3.2 Emission maps at $1.2\,\rm mm$

Fig. 4 plots the projected $1.2\,\rm mm$ emission maps of our models $1$-$6$. The emission intensity is normalised, the black color representing the faintest pixels and the white one the brightest ones, respectively. The white bar in the bottom right corner of each figure indicates the physical size of each images. The top series of panels represent the young disc at times $20.1\,\rm yr$ (Fig. 4a), $23.0\,\rm yr$ (Fig. 4b) and $26.4\,\rm yr$ (Fig. 4c) after the beginning of the simulation, respectively. These snapshots correspond to time instances before the MYSO experiences its first accretion-driven burst. The disc model 1 in Fig. 4a is a circular and stable disc that does no show any noticeable sign of either fragmentation by gravitational instability or substructures in it. This shape persists up to at least $23.0\,\rm yr$ (model 2, Fig. 4b). During this time interval, the protostellar mass evolves from $8.6,\rm M_{\odot}$ to $10.7\,\rm M_{\odot}$, but its circumstellar medium fails in developing a clear spiral structure, although a faint trailing structure begins to appear at radii $\geq 200$-$600\,\rm au$. The brightest region is the inner $50\,\rm au$ of the disc. Our model 3 (Fig. 4b) starts exhibiting a two-arms pattern extending up to radii $\sim 400\,\rm au$ from the MYSO. Moreover, a dense gaseous clump formed in the spiral arm inward-migrates and fast-moves onto the protostar. The clump keeps on contracting during its migration and its core increases in density and temperature, constituting the maximal-emitting region of the disc. The successive, rapid accretion of several migrating clumps from different parent arms, the merging of clumps into inhomogeneous gaseous structures or even the migration of clumps that separate an inner portion of spiral arms into two segments make the accretion pattern more and more complex and strengthen the accretion variability (Fig. 1).

The bottom series of panels of Fig. 4 shows emission maps of our disc models 4 to 6, corresponding to the young disc at times $30.7\,\rm yr$ (Fig. 4c), $31.1\,\rm yr$ (Fig. 4e) and $32.1\,\rm yr$ (Fig. 4f), respectively. The disc model 4 (Fig. 4d) has the typical morphology of the circumstellar medium of a protostellar disc undergoing efficient gravitational instability. In addition to clear enrolled spiral arms starting close to the star and fading away at $\geq 600$-$800\,\rm au$, several gaseous clumps have formed inside the spiral arms, at distances $\geq 200$-$400\,\rm au$ from the protostar. The spiral arms and clumps appear very prominently in the images, with a higher integrated surface brightness than the homogenous parts of the disk. The disc model 5 (Fig. 4e) is more complex than the rather standard ring-like clump-hosting structure of model 4. Its strongly fragmented disc experiences a complex dynamics in which multiple Toomre-instable gaseous clumps (Meyer et al., 2018) exert torques with each other, with the neighbouring spiral arms and with the central MYSO of $16.7\,\rm M_{\odot}$, respectively. Finally, the disc model 6 (Fig. 4f) exhibits an extended, strongly fragmented disc with a filamentary spiral arm hosting several bright gaseous clumps. This morphological configuration of the disc substructures is typical of the time following an accretion-driven burst, when the tail of a migrating clump is swung away and pushes the circumstellar material to the larger radii as described in Meyer et al. (2017). The subsequently forming dense ring is the location of the formation of new clumps, whose integrated surface brightness can be similar to that of the inner disc.

3.3 Observability of evolving discs around growing protostars

Figs. 5 and 6 show emission maps of our disc models as seen by alma if located at a distance of $1\,\rm kpc$ from the observer and monitored with antenna configuration 8 and 10, respectively. Both figures show $1.2\,\rm mm$ simulated images of the selected snapshots at times $20.1\,\rm kyr$ (a), $23.0\,\rm kyr$ (b), $26.4\,\rm kyr$ (c), $30.7\,\rm kyr$ (d), $31.1\,\rm kyr$ (e) and $32.1\,\rm kyr$ (f), respectively. The images assume a face-on viewing geometry. On each panel the beam size is given in the bottom left panel. The map scale is in units of $\rm arcsec$ and the surface brightness intensity is in $\rm mJy\,\rm beam^{-1}$. The synthetic images of stable disc model 1 exhibit a bright circular shape for both antenna configuration 8 (Fig. 5a) and the most extended configuration 10 (Fig. 6a). The same is the case for our disc model 2 ($23.0\,\rm kyr$), as the current beam resolution of alma is not able to spatially resolve the early forming overdensities in its principal spiral arm (Fig. 5b and Fig. 6b). The $1.2\,\rm mm$ signature of the effects of gravitational instability begins to appear at time $26.4\,\rm kyr$ with the antenna configuration 10 (Fig. 6c), when the protostar has reached $12.1\,\rm M_{\odot}$ and a series of twin spiral arms in which a detached migrating clump falls onto the protostar. The images at time $30.7\,\rm kyr$ clearly reveal the fragmenting characteristic of our disc, i.e. a self-gravitating spiral arm that includes a few bright gaseous clumps (Fig. 5d) potentially on the way to low-mass star formation (Meyer et al., 2018). The accretion disc and its substructures can be resolved with the alma antenna configuration 8. A more extended antenna configuration (configuration 10) yields similar results, see Fig. 6d. At time $\geq 31.1\,\rm kyr$ the synthetic observations of models 5 and 6 reveal a nascent multiple, hierarchical massive system with clumps and filamentary structures organised around a hot, young massive star of mass $16.7\,\rm M_{\odot}$ (Figs. 5e,f) and $20.0\,\rm M_{\odot}$ (Figs. 6e,f), respectively.

Shown in Fig. 7 are cross-sections extracted from the alma cycle 6 maps with antenna configuration 8 (Fig. 5) and antenna configuration 10 (Fig. 6), respectively. The figures compare the emission intensity at times $20.1\,\rm kyr$ (thin solid blue line), $23.0\,\rm kyr$ (solid blue line), $26.4\,\rm kyr$ (thick solid blue line), $30.7\,\rm kyr$ (thin dashed red line), $31.1\,\rm kyr$ (dashed red line) and $32.1\,\rm kyr$ (thick dashed red line). These profiles assume a face-on viewing geometry ($\phi=0\mathrm{\SIUnitSymbolDegree}$). The source distance is $1\,\rm kpc$ and the cuts is taken along the vertical direction through the center of the synthetic images. The shape of the cross-sections (solid blue lines) corresponding to Fig. 5 further illustrates that the surface brightness of our modelled accretion disc does not exhibit signs of substructured circumstellar disc at times $\leq 23.0\,\rm kyr$, when the protostar is $\leq 10.7\,\rm M_{\odot}$, although the intensity cut extracted from the image at $26.4\,\rm kyr$ (thick solid blue line) has two additional bright emission peaks corresponding to the two spiral arms.

The slightly negative fluxes in some places in the interarm region of the discs result from the noise of the substracted background thermal emission arising from the incomplete u-v coverage in the simulated observations.

The cuts taken in the models 4-6 clearly highlight the presence of a dense ring-like spiral arm that surrounds the MYSO and which produces intensity peaks of $\approx 75\,\rm mJy/beam$ above the background of $\approx 10\,\rm mJy/beam$.

The emission cuts in Fig. 7b have a higher spatial resolution but a lower sensitivity than those of Fig. 7a, i.e. the intensity differences between the close stellar surroundings and the arms/clumps in the disc are less pronounced (see also Fig. 14).

For example, the northern ring of model 6 ($32.1\,\rm kyr$) peaks at $\approx 50$ and $\approx 30\,\rm mJy/beam$, which represents about $25\%$ and $60\%$ of the disc central emission, respectively. The images in Fig. 5a-f are complementary to that of Fig. 6a-f in the observation and analysis of the discs, and observations of the circumstellar medium of a high-mass star should be performed with several antenna configurations of alma for a given waveband.

3.4 Effect of the inclination angle

Figs. 8 and 9 show synthetic alma cycle 6 ($1.2\,\rm mm$) images of our disc models. In the figure, each row represents an inclination, and each column a time, and the images are calculated for the alma C43-8 configuration, assuming a distance to the source of $1\,\rm kpc$. The beam size is given in the lower left corner of each panel. Again, the map scale is in units of $\rm arcsec$ and the surface brightness intensity is in $\rm mJy\,\rm beam^{-1}$. For $\phi=45\mathrm{\SIUnitSymbolDegree}$, the emission maps do not exhibit important changes in the disc shapes, except that their projected size is squeezed by a factor $\cos(\pi/4)$ along the direction normal of the disc plane (the vertical direction of the maps). This is consistent with the studies on discs from low-mass stars of Dong et al. (2016) which demonstrated that moderate inclination angles with respect to the plane of sky have little incidence in the disc emission maps compare to that with $\phi=0\mathrm{\SIUnitSymbolDegree}$. In the case of edge-on accretion discs, the antenna configuration 8 does not allow the observer to clearly distinguish clumpy structures in the accretion discs, even if they can be inferred from the meridional emission profile, e.g. at time $32.1\,\rm kyr$ (Fig. 8l). Similar trends in the morphological signatures of the accretion discs happen with a higher resolution of the beam, when the antenna configuration is more extended (Fig. 9).

Fig. 10 compares cross sections taken from our simulated alma cycle 6 ($1.2\,\rm mm$) images, with antenna configuration 8 (compact, coarser resolution) and 10 (extended, higher resolution). The source distance is $1\,\rm kpc$ and all cuts are taken along the horizontal direction through the location of the star. The first panel Fig. 10a illustrates the main difference between face-on and edge-viewed accretion discs during the early phase, at times $\leq 26.4\,\rm kyr$. While the face-on disc ($\phi=0\mathrm{\SIUnitSymbolDegree}$) has first an emission profile that does not exhibit clear substructures, its edge-view ($\phi=90\mathrm{\SIUnitSymbolDegree}$) indicates the presence of its large, developing spiral arm. At later time and by projection, it gives the disc a symmetric morphology but not a uniform surface brightness with respect to its axis of rotation.

Similar conclusions are drawn with the antenna configuration 10 (Fig. 10c). The disc projected intensity deceases as a function of the inclination, as other circumstellar structures around more evolved massive stars do (Meyer et al., 2016). The time evolution of the surface brightnesses of the edge-viewed accretion discs is as follows: the intensity peaks on the projection plane at the location of the star. The emission signature gradually increases in size as the disc grows, by accumulating molecular gas from the surrounding infalling envelope, up to reaching a radius of $\simeq 800$-$1000\,\rm au$ (thick red curves of Fig. 10b,d). As a result of the antenna configuration 10 and smaller beam size, the meridional cuts through the disc reveal more details in the projected disc emission. Interestingly, some intensity peaks such as at time $30.7\,\rm kyr$ (thin black dotted dashed line of Fig. 10d) reveal the presence of a couple of aligned clumps, while the other snapshots do not. Polarization measurements in the near-infrared waveband may be useful to further investigate this phenomenon.

3.5 Effect of the distance of the source

Although the Earth’s closest high-mass star-forming regions are located at distances $\leq 1\,\rm kpc$ (Reid et al., 2009), numerous giant molecular regions with massive protostars are more distant. We check the observability of the substructures in our accretion disc models assuming a distance of $\leq 2\,\rm kpc$. The radiative transfer methods and the computations of synthetic images are identical as above described for sources located at $1\,\rm kpc$, the only difference being the exposure time that we derive to be $30\,\rm min$ (Section 4.2). Fig. 11 shows synthetic alma cycle 6 ($1.2\,\rm mm$) images of our disc models.

In the figure, each row represents an inclination, and each column a time, and the images are calculated for the alma C43-10 configuration, assuming a distance to the source of $2\,\rm kpc$.

The beam size is given in the lower left corner of each panel. Again, the map scale is in units of $\rm arcsec$ and the surface brightness intensity is in $\rm mJy\,\rm beam^{-1}$. In accretion discs located at $2\,\rm kpc$ rather than at $1\,\rm kpc$, it is more difficult to observe clearly detached clumps for young discs such at times $\leq 26.4\,\rm kyr$ whatever their viewing geometry is (Fig. 11a,e,i). Synthetic images of later time instances ($>26.4\,\rm kyr$) lead to similar conclusions regarding the possible observability of clumps and spiral arms as those derived on the basis of the synthetic images computed assuming a distance of $1\,\rm kpc$. In the case of edge-on discs, inferring clumps is somewhat less evident (Fig. 11i-k), even though the far West clump developing at time $32.1\,\rm kyr$ can be distinguished from the rest of the circumstellar structure in Fig. 11l. We conclude that the outcomes of gravitational instabilities in accretion discs around MYSOs can be detected with alma, with observations performed during the Cycle 7 observational run and using its most extended available antenna configuration.

4 Discussion

This last section is devoted to the discussion of our synthetic observations of fragmenting circumstellar discs around MYSOs. We first review the limitations and caveats of our methods, investigate how the effects of the integration time of the interferometric observations and of the alma facility beam resolution can affect the detection of substructures in the disc. We briefly compare our results with precedent studies and we discuss future prospects for observing gravitational instability and disc fragmentation around MYSOs. Finally, we discuss our results in the context of high-mass protostellar disc-hosting candidates.

4.1 Limitations of the model

Our models of synthetic images are affected by two kinds of caveats: on the one hand the numerical hydrodynamics simulations of our forming accretion discs are intrinsically subject to limitations justifying future improvements, on the other hand the radiation transfer methods that we use in this study might, in their turn, be improved as well. The principal limitation of our models comes from the restrictions on the hydrodynamical timestep, controlled by the standard Courant-Friedrichs-Levy parameter (Mignone et al., 2007, 2012). It decreases proportionally to the smallest grid size of the grid mapping the computational domain, and, therefore gets smaller as the spatial resolution augments. The sink cell radius also influences the speed of the time-marching algorithm in our simulations. It should be as small as possible since we study the accretion of inward-migrating disc fragments by the protostellar surroundings. In this study, we adopt the smallest value ($r_{\rm in}=10\,\rm au$) where the simulations are still computationally affordable, which is still smaller than in recent works on disc fragmentation. Nevertheless, the most critical point of the disc simulations in spherical coordinate systems remains the effect of the inertia of the central star to the disc dynamics, which is neglected when the high-mass protostar is fixed at the origin of the computational domain, see also discussion in Meyer et al. (2019). The stellar motion in response to the gravitational force of the disc has been shown not only to reduce the strength of gravitational instability (but not shutting it down completely) in the low-mass regime of star formation but also to introduce unpleasant boundary effects (Regály & Vorobyov, 2017), as it displaces the barycenter of the star-disc system from the origin of the domain. This is in accordance with the analytic study of Adams et al. (1989) that pointed out the role of stellar inertia in the development of asymmetries in discs. The study of Meyer et al. (2019) demonstrated that the disc inertial force does not entirely stabilise discs or shuts off the burst phases of the protostellar growth, in accordance with results of Regály & Vorobyov (2017). The assumptions in our numerical simulations of the disc structures are presented and discussed in detail in Meyer et al. (2018) and Meyer et al. (2019).

The dust temperature is a key parameter governing the brightness of the accretion discs. Fig. 12 compares simulated ALMA cycle 6 (1.2 mm) images with an antenna configuration 10 of a stable (left) and a fragmented (right) disc model, where the radiation transfer calculations are performed with different methods to compute the dust temperature fields: either by Monte-Carlo in the radmc-3d code (top panels) or using the dust temperatures from the hydrodynamics simulations (bottom panels). The images are shown for face-on viewing angle ($\phi=0\mathrm{\SIUnitSymbolDegree}$). The beam size is indicated in the lower left corner of each panel and the source distance is $1\,\rm kpc$. It appears that the synthetic alma images of the stable model at time $20.1\,\rm kyr$ do not exhibit strong differences if modelled on the basis of either of a Monte-Carlo method or using the temperature fields imported from the hydrodynamical simulations (left panels of Fig. 12). Nevertheless, differences are observed with respect to the observability of the clumps of fragmented discs (right panels of Fig. 12). Such changes in the simulated images are due to the fact that the hydrodynamical temperature field results in a self-consistent simulation, including the treatment of the mechanical $pdV$ work from gravity forces in the disc, i.e. the internal thermodynamics of the clumps are taken into account. The cores of these clumps experience an internal gravitational collapse (Meyer et al., 2018), while a Monte-Carlo calculation of the temperature fields is less realistic as it simply ray-traces photon packages through the clumps without accounting for their thermal properties. This results in brighter $1.2\,\rm mm$ clumps in the spiral arms of the fragmented discs (bottom right panel of Fig. 12), and, consequently, in a higher brightness and better observability at the considered sub-millimetre wavelength. This effect was also reported in Dong et al. (2016) in the context of fragmented disks around low-mass protostars.

Fig. 13 compares several of the cross-sections taken through our simulated ALMA cycle 6 ($1.2\,\rm mm$) antenna configuration 10 images displayed in Fig. 12. Two disc models are compared: a first one corresponding to a time $20.1\,\rm kyr$ when the disc is stable (blue lines) and another one at time $30.7\,\rm kyr$ when the disc fragments (red), in which the dust temperature is estimated by Monte-Carlo simulations (thin solid lines) and by thermal equilibrium (thick dashed lines).

It further highlights the differences between the two methods for the determination of the disc dust temperature. One can see that the simulated images performed with two consecutive Monte-Carlo ray-tracings (RADMC-3D) are slightly fainter in intensity than that carried out on the basis of PLUTO temperatures, e.g. the central part of the RADMC-3D image peaks at $\approx 35\,\rm mJy/beam$ while the PLUTO peaks at $\approx 50\,\rm mJy/beam$, respectively. As another example, the peaks intensities in the clumps calculated with RADMC-3D and PLUTO temperatures are $\approx 10$ and $\approx 14\,\rm mJy/beam$, respectively, which is slightly more noticeable than at the location of the central high-mass protostar. Our simulations with radmc-3d tends to underestimate the surface brightness in the synthetic images. As a conclusion, we stress the need for the use of self-consistently determined temperatures for the modelling of emission maps of accretion discs of massive protostars.

4.2 Effect of integration time and angular resolution

Fig. 14 compares the effects of the antenna configuration and of the integration time on the simulated images of our fragmented accretion disc at time $30.7\,\rm kyr$. The top panels are simulated with the most compact antenna configuration (cycle 6.6) and lower panels are simulated with the most extended configuration (cycle 6.10), respectively. The integration time is $600\,\rm s$ (left) and $1800\,\rm s$ (right). All images are shown assuming a face-on viewing angle of the discs. The beam size is indicated in the lower left corner of each panel and the source is located at $1\,\rm kpc$. The series of tests highlight that the more extended the antenna configuration, i.e. the smaller the beam scanning a region of the sky, the better visible are the substructures inside the protostellar disc. If the images performed with $600\,\rm s$ of integration time and an antenna configuration 6 do not reveal clumps in the disc, those with antenna configuration $\geq 7$ do. The latter show both the distinct features of two clumps which are brighter than the background, and it equivalently reveals the main spiral arm of the disc. Note that the best available beam resolution (bottom panels) even shows the asymmetric character of some clumps. The two series of panels comparing the alma integration time of $600\,\rm s$ (right) and $1800\,\rm s$ (left) indicate that the more extended the antenna configuration, the faster the convergence of the solution in terms of exposure time. Indeed, the more compact antenna configuration 6 (top series of panel) shows that the image modeled with $1800\,\rm s$ reveals a clump above the protostar that is not visible with only $600\,\rm s$ integration time, while the other images with the antenna configuration $\geq 7$ do not exhibit such differences, i.e. the beam resolution is more important than the integration time. Similarly, we derive a reasonable exposure time to observe sources located at $2\,\rm kpc$ to be $30\,\rm min$.

4.3 Comparison with other studies

A couple of previous works have performed synthetic imaging of both analytic and self-consistent hydrodynamical simulations of forming massive stars at the disc scale ($\leq 1000\,\rm au$), see Krumholz et al. (2007), Harries et al. (2017), Meyer et al. (2018) and Jankovic et al. (2019). First, the paper of Krumholz et al. (2007) investigated synthetic alma and noema images of disc structures forming from collapsing turbulent pre-stellar cores located at $\leq 2\,\rm kpc$, in which secondary star formation is tracked with the used of Lagrangian sink particles whose caveats are discussed in Meyer et al. (2018). Our models therefore provide more accurate density and temperature disc structures, as the gaseous clumps and secondary low-mass stars are self-consistently produced and their internal thermodynamics is intrinsically taken into account. Secondly, Harries et al. (2017) carried out state-of-art gravito-radiation-chemo-hydrodynamical simulations and checked the disc spiral arm observability with simulated $1\,\rm mm$ alma interferometric observations. This study includes the most accurate treatment of the radiative protostellar feedback to-date, however, its disc models neither strongly fragment nor produce the formation of clumps because of a lack of spatial resolution.

In the earlier studies of this series devoted to the burst mode of massive star formation, Meyer et al. (2018) showed that spiral arms can be observed in continuum and molecular emission for a fragmenting disc located at a distance and with a inclination corresponding to the disk around the MYSO AFGL-4176 (Johnston et al., 2015). Last, the work of Jankovic et al. (2019) predicts the observability of spiral arms and clumps in discs around massive stars by computing images of analytic disc structures. They predict their observability up to a distance $\leq 5\,\rm kpc$ for arbitrary inclinations by direct imaging discs with $\phi\leq 50\mathrm{\SIUnitSymbolDegree}$ and/or by their kinematic signature. Our results show that clumps can, under some circumstances, be directly imaged or at least inferred for edge-viewed discs. These studies highlight the possible observation of circumstellar discs around young massive stars. Our work allows for the first time to obtain self-consistent solution from the hydrodynamical, stellar evolution and radiative transfer point of view, which permits qualitative predictions relative to the emission properties of disc substructures.

4.4 Prospect for observing disc fragmentation around MYSOs: the cases of S255IR-NIRS3 and NGC 6334I-MM1

In this sub-section, we discuss the prospects for the potential observability of signatures of gravitational instabilities in the surroundings of outbursting sources. We focus on the two high-mass protostars for which potentially accretion-driven bursts have been observed, namely S255IR-NIRS3 (Caratti o Garatti et al., 2017) and NGC 6334I-MM1 (Hunter et al., 2017). Note that a third burst from a MYSO, V723 Car, has been reported in Tapia et al. (2015), however, its discovery in archival data did not permit further monitoring of it. S255IR-NIRS3 is the very first massive protostar observed experiencing a disc-mediated burst and it is therefore a natural target to investigate the possible presence of other clumps in its disc. NGC 6334I-MM1 is the second known bursting MYSO and it exhibited a successive quadruple temporal variation in the dust continuum emission which suggests that this eruptive event is stronger in intensity than that of S255IR-NIRS3. This makes NGC 6334I-MM1 the second most obvious candidate to chase for substructures in its close surroundings. A summary of the recent research concerning these two MYSOs is given in Section 5.4 of Meyer et al. (2019). The predictions for their circumstellar medium are performed following the comparison method utilised for the proto-O-type star AFGL 4176 in Meyer et al. (2018). We select a simulation snapshot of our hydrodynamical model at a time corresponding to the estimated masses of the MYSOs and run radiative transfer calculations as above described, with radmc3d and casa. The dust continuum models are run using the values for $L_{\star}$, $R_{\star}$ and $T_{\rm eff}$ either constrained by observations or derived from our stellar evolution calculation, and we account for the position of the MYSO on the sky, its distance to the observer and, when possible, the inclination angle of the disc with respect to the plane of sky.

We adopt the stellar parameters of S255IR-NIRS3 and NGC 6334I-MM1 summarised in the studies of Caratti o Garatti et al. (2017), Eislöffel & Davis (1995), Bøgelund et al. (2018) and Hunter et al. (2017). For S255IR-NIRS3, we consider that it has recovered the pre-burst quiescent value of $L_{\star}$, although our stellar evolution calculation predict a rather slow post-burst luminosity decay lasting $\simeq 10^{2}\,\rm yr$, see our Fig. 2 and Meyer et al. (2019). Our choice to adopt protostellar parameters fainter than that of Run-1-hr is justified because we have classified the flare of S255IR-NIRS3 to be a modest 2-mag burst in Meyer et al. (2019) which can not be responsible for long-lasting spectral excursions in the Hertzsprung-Russell diagram as experienced by the MYSO of Run-1-hr, see also Meyer et al. (2019). Concerning NGC 6334I-MM1 which experienced a strong burst, one has to keep in mind than the values of $R_{\star}$ in Forgan et al. (2016) do not account for the burst-induced protostellar bloating. However, the values of $T_{\rm eff}$ are fully consistent with our stellar evolution calculation. We correct this discrepancy by assuming a larger value of the protostellar radius, taken from Run-1-hr, which in its turn also affects $L_{\star}$. Hence, the only unknown is the inclination angle $\phi$ of the possible disc of NGC 6334I-MM1 with respect to the plane of sky. The various parameters of these two MYSOs are summarised in Tab. 2 and our results are shown in Figs. 15 and 16.

The Cycle 7 band 6 synthetic observation of S255IR-NIRS3 is plotted in Fig. 15. The disk position angle is chosen to match the direction of the two lobes inside which material is expelled during the burst in Caratti o Garatti et al. (2017). It is interesting to see that a few clumps in the disc can still be observed despite of the almost edge-on view of the disc around S255IR-NIRS3. The Cycle 7 band 6 synthetic observation of NGC 6334I-MM1 is shown in Fig. 16. As the inclination of its disc is not known, we assumed several angles ($\phi=0\mathrm{\SIUnitSymbolDegree}$, $30\mathrm{\SIUnitSymbolDegree}$ and $90\mathrm{\SIUnitSymbolDegree}$). We conclude for all the considered angles $\phi$ that substructures can easily be observed in the disc and its surroundings. Hence, these two targets are good candidates for further exploration of substructures in their circumstellar discs. Our results also permit to speculate that NGC 6334I-MM1 might be recovering its quiescent, luminous ground state after a spectral excursion in the red part of the Hertzsprung-Russell. Consequently, this object is also a candidate to test our prediction regarding the intermittency of the H ii  region surrounding outbursting MYSOs that was predicted in Meyer et al. (2019).

One should keep in mind that the two MYSOs that we selected to study the observability of their (possibly substructured) accretion discs are the two MYSOs which had a known burst. This was intentional, as the burst mode picture of star formation predicts a direct correspondance between the occurrence and properties of the burst with the effects of gravitational instability in those discs (Meyer et al., 2017). However, all MYSOs may have an accretion disc and the absence of observed bursts in most of them can simply be interpreted as them being in a quiescent phase of accretion, not as the absence of discs or even not as the absence of fragments in their discs. Indeed, the numerical models for bursting massive protostars of Meyer et al. (2019) indicate that MYSOs spend $<2\,\%$ of their early $60\,\rm kyr$ being in the eruptive state. Consequently, the other young massive stars with candidate discs might also be worth observing with future high-resolution alma campaigns. Possible targets are listed in the observational studies reporting the Red MSW Source survey of Galactic hot massive protostars (Lumsden et al., 2013). The closest ones, located at $\approx\,1.4\,\rm pc$ are G076.3829-00.6210 (S106) and G080.9383-00.1268 (IRAS 20376+4109) while a few others are at a distance of $\approx\,2.0\,\rm pc$. For instance, G015.0357-00.6795 (M17), the W3IRS objects, G192.5843-00.0417 (G192.58) might deserve particular attention for the observer observational astronomical concentrating on multi-epoch VLBI (very long baseline interferometry) observations of maser emission.

Last, we should also mention that there is a growing suspicion of correlations between the maser emission from MYSOs and their recent accretion activity, as maser flares originating from the surroundings of massive protostars continue to be monitored (Burns et al., 2016; Burns et al., 2017; MacLeod et al., 2018). This promising research avenue might investigate, amongst others, if maser variability is an outcome of episodic accretion. Such results would be of great interest not only in the understanding of massive star formation mechanisms but also a huge leap forward in the validity of the burst mode picture of star formation.

5 Conclusions

This work explores the observability of a fragmenting accretion disc surrounding massive young stellar objects (MYSOs), as modelled in Meyer et al. (2017, 2018), using the alma interferometer. Such discs are subject to efficient gravitational instabilities, generating the formation of dense spiral arms and gaseous clumps, which episodically inward-migrate to the inner disc region. Some of them further contract up to become the precursors of close/spectroscopic binary companions to the central protostar (Meyer et al., 2018) and/or (simultaneously) provoke accretion driven outbursts of various intensities (Meyer et al., 2019). Representative simulation snapshots exhibiting the characteristic formation and evolutionary stages of the disc are selected from model run-1-hr of Meyer et al. (2018) which results from the gravitational collapse of a $100\,\rm M_{\odot}$ molecular pre-stellar core. We performed stellar evolution and radiative transfer calculations to produce alma at band 6 ($1.2\,\rm mm$) emission maps of the discs and obtain synthetic images. These images are then used to investigate the observability of spiral arms and clumpy fragments. In addition to the disc evolution, we explore how its inclination angle with respect to the plane of sky affects the $1.2\,\rm mm$ projected emission signatures. We also look at the effects of various alma beam resolutions and/or antenna configurations to predict the best facility configuration and the ages of accretion discs from MYSOs that are the easiest ones to observe.

Our results demonstrate that not only the spiral arms developing in the disc by gravitational instability, but also the gaseous clumps can be observed, as their internal density and temperature make them much brighter than the disc itself. For sources located at a distance of $1.0\,\rm kpc$, $10\,\rm min$ exposition time with a $0.015^{\prime\prime}$ beam resolution permits to clearly distinguish those nascent stars and reveal the multiplicity in the environment of young massive protostars. We predict that in the millimeter waveband (band 6), an antenna configuration 8 (corresponding to a beam resolution of $0.029^{\prime\prime}$) is necessary to reach a beam resolution required to identify the disc substructures with a rather short $10\,\rm min$ exposure time. However, only antenna configuration 10 (corresponding to a beam resolution of $0.015^{\prime\prime}$) might allow us to distinguish the finest disc substructures. Interestingly, our alma images also revealed clumps in inclined discs and they can be inferred for edge-on accretion discs, as long as these fragments develop at radii $\geq 600\,\rm au$ from the protostar, i.e. they are not fully screened by the irradiation of the central massive protostar. In our models, this corresponds to discs older than $\approx 26.7\,\rm kyr$ after the onset of the gravitational collapse, when the MYSO has already reached a mass of $\approx 16\,\rm M_{\odot}$. Hence, good candidates for hunting disc substructures are MYSOs of $\simeq 20\,\rm M_{\odot}$. Our conclusions regarding the observability of substructures in our disc models are valid up to distances of $2\,\rm kpc$, using a $0.015^{\prime\prime}$ beam resolution with the alma Cycle 7 antenna configuration 10 and at least $30\,\rm min$ exposure time.

The detection probability of the structures increases as the disc grows, i.e. accretion discs around protostars close to the ZAMS are more extended, and host more clumps in the outer region of the disc that is not screened by the MYSO irradiation. Infalling clumps in the inner few $~{}10-100\,\rm au$ of the disc might not by observable and the probability to directly associate a clump with a burst is rather small. Last, our study emphasises the necessity to perform simulated images of discs around massive stars by accounting for the gaseous clumps thermodynamics by using temperatures self-consistently determined during the hydrodynamics simulations instead of dust temperature fields obtained using steady-state radiative transfer models which neglect $pdV$ work. In addition to demonstrating the possible direct imaging of nascent multiple massive protostellar systems with alma, we also highlight that gaseous clumps and spiral overdensities can clearly be observed/inferred at distances corresponding to the bursting massive protostars S255IR-NIRS3 (Caratti o Garatti et al., 2017) and NGC 6334I-MM1 (Hunter et al., 2017). Those objects might deserve particular future attention from the point of view of high-resolution observations (Zinchenko et al., 2015). Similarly, the Earth’s closest sources, at $\approx\,1.4-2.0\,\rm pc$, of the Red MSW Source survey of Galactic hot massive protostars (Lumsden et al., 2013) also constitute valuable potential targets. As the models for fragmented discs of Meyer et al. (2017, 2018) predict a pre-main-sequence outburst of massive protostars to be the direct outcome of rapidly inward-migrating clumps, we suggest observers to consider their respective circumstellar discs as priority targets to search for the substructures which can be signatures of fragmentation by gravitational instability in their accretion disc.

Acknowledgements

The authors thank the referee for her/his valuable comments which improved the quality of the manuscript. D. M.-A. Meyer acknowledges B. Magnelli from the German ALMA Regional Center (ARC) Node at the Argelander-Institut für Astronomie of the University of Bonn for his help, and C. Dullemond for his advices in using the radmc-3d code. D. M.-A. Meyer also thanks Dr. Z. Regály from Konkoly Observatory for his advices in using the casa software and Dr. K. G. Johnston for kindly sharing her explanations with the matplotlib python package. S. Kraus and A. Kreplin acknowledge support from an ERC Starting Grant (Grant Agreement No. 639889), a STFC Rutherford Fellowship (ST/J004030/1) and an STFC Rutherford Grant (ST/K003445/1). E. I. Vorobyov acknowledges support from the Russian Science Foundation grant 18-12-00193. This work was sponsored by the Swiss National Science Foundation (project number 200020-172505).