Radiation-hydrodynamical simulations of massive star formation using Monte Carlo radiative transfer II. The formation of a 25 solar-mass star

TL;DR

This study uses Monte Carlo radiation hydrodynamics simulations to model the formation of a 25-solar-mass star, revealing accretion processes, disc structure, and observable features consistent with real observations.

Contribution

It presents a detailed RHD simulation of massive star formation with synthetic observations, advancing understanding of accretion and disc dynamics in such systems.

Findings

Star reaches 25 solar masses after 35,000 years

Disc mass is 7 solar masses with near-Keplerian velocity

Synthetic IR images align with observations, but some discrepancies exist

Abstract

We present a numerical simulation of the formation of a massive star using Monte- Carlo-based radiation hydrodynamics (RHD). The star forms via stochastic disc accretion and produces fast, radiation-driven bipolar cavities. We find that the evolution of the infall rate (considered to be the mass flux across a 1500 au spherical boundary), and the accretion rate onto the protostar, are broadly consistent with observational constraints. After 35kyr the star has a mass of 25 solar masses and is surrounded by a disc of mass 7 solar masses and 1500 au radius, and we find that the velocity field of the disc is close to Keplerian. Once again these results are consistent with those from recent high-resolution studies of discs around forming massive stars. Synthetic imaging of the RHD model shows good agreement with observations in the near- and far-IR, but may be in conflict with observations that suggests that MYSOs are typically circularly symmetric on the sky at 24.5 microns. Molecular line simulations of a CH3CN transition compare well with observations in terms of surface brightness and line width, and indicate that it should be possible to reliably extract the protostellar mass from such observations.

Tables

Table 1. Table 1: Observed properties of massive protostellar disc candidates along with the model disc parameters from this study. The columns give the object designation, the estimated mass of the disc, its estimated radius, the approximate disc-to-star mass ratio, and the reference.

Object	Mass ( M⊙)	Radius (au)	Mdisc/M∗	Reference
NGC 6334 I(N), SMA 1b	10–30	400	$\sim 0.5$	Hunter et al. (2014)
IRAS 16547–4247	6	500	$\sim 0.5$	Zapata et al. (2015)
G11.92-0.61 MM1	2–3	1200	$\sim 0.05$	Ilee et al. (2016)
AFGL 4176	12	2000	$\sim 0.5$	Johnston et al. (2015)
IRAS 20126+4104	1.5	858	$\sim 0.13$	Chen et al. (2016)
Our model	7	1500	0.3	This work

Equations

$$X_{\rm CH_{3}CN}=\left\{\begin{array}[]{ll}10^{-8}&T>100\,{\rm K}\\ 5\times 10^{-9}&90\,{\rm K}<T<100\,{\rm K}\\ 10^{-10}&T<90\,{\rm K}\end{array}\right.$$

Full text

††pagerange: Radiation-hydrodynamical simulations of massive star formation using Monte Carlo radiative transfer II. The formation of a 25 solar-mass star–Acknowledgments††pubyear: 2016

1 Introduction

Although the paradigm for low-mass star formation has achieved a broad consensus, the mechanism by which high mass stars form is rather less clear (e.g. Tan et al. 2014). Observationally this is because the rapidity and rarity of massive star formation means that the process occurs while the stars are deeply embedded and that we observe it at great distances. However rapid progress is now being made, with the advent of high resolution observations of massive young stellar objects (MYSOs) using ALMA, that encode geometrical, thermal, chemical and dynamical information on the protostars, discs, and envelopes–for a recent review see Beltrán & de Wit (2016).

From a theoretical perspective massive star formation is challenging because theorists must deal not only with the hydrodynamics of self-gravitating collapse, but also with the enormous radiation field produced by the massive protostar, which not only heats and ionises the gas, but also results in forces that act to resist the collapse, leading to outflows and disc winds. Progress in numerical simulations has followed a pattern of increasing complexity in the treatment of these radiative-feedback processes. The first calculations relied the diffusion approximation, which is easy to implement and computationally cheap and captures the physics of the radiation field well in the optically thick limit. The diffusion approximation breaks down as the radiation-field becomes free-streaming, and a flux-limiter is used in optically thin regions to overcome this drawback. Flux-limited diffusion (FLD) has been used successfully to model high-mass star formation (Yorke & Sonnhalter, 2002; Krumholz et al., 2009).

One potential pitfall of using grey FLD to study high-mass star formation is the enormous disparity between the shape of the radiation field of the protostellar sources (often characteristic of temperatures exceeding 10 kK) and that of the thermal dust emission ($<1500$ K). The fact that dust absorption and scattering cross sections are much larger in the blue means that grey FLD may significantly underestimate the radiation forces from the direct (unprocessed) protostellar radiation field. This drawback is avoided by the so-called hybrid RT scheme, in which a frequency-dependent ray-tracing integral of the RT equation from the protostellar source is used to determine the intensity of the radiation field until it becomes optically thick, and the radiation field at this point is used as a boundary condition for the diffusion approximation in the rest of the computational domain. This method was first used in 1-D simulations of massive star formation by Wolfire & Cassinelli (1986, 1987) and was extended to higher dimensions by Kuiper et al. (2010). Although more computationally expensive than FLD, the method has the advantage that it captures the physics of the first absorption of the protostellar radiation field by the dust, as well as any shadowing effects close to the source. Tests of the method show that it does a remarkably good job of determining the dust temperature when compared to dedicated RT codes (Kuiper & Klessen, 2013). Many contemporary studies now employ this hybrid method (Kuiper et al., 2012; Klassen et al., 2016; Rosen et al., 2016).

Of course one would really like to conduct the RT step of the RHD simulation using a method that treats both the dust and gas opacities and emissivities in their full polychromatic glory. Dedicated RT codes such as cloudy (Ferland et al., 2013) manage this, as do a large number of Monte Carlo (MC) based codes (Bjorkman & Wood, 2001; Harries et al., 2004; Wood et al., 2004; Ercolano et al., 2005; Pinte et al., 2006; Robitaille, 2011), but these are usually perceived to be too slow to be useful as the radiation step in an RHD code. However the MC method has one distinct advantage: photon packets are independent. This means that MC loops are embarrassingly parallelisable, and one can use computer science to exploit the full power of modern distributed supercomputers to bring down the overhead of a MC RT step to something approaching that of the hydrodynamics step itself (Haworth & Harries, 2012). Furthermore once the simulation is complete the temperature of the dust and gas, and the photoionisation state of the gas may be used used to post-process the simulation and create synthetic observables, allowing a direct comparison with observations (Haworth et al., 2012).

In Harries (2015) we presented a new MC method for including radiation pressure in RHD simulations that incorporates a level of microphysical detail that is significantly greater than that of flux-limited diffusion or hybrid techniques. We have shown that the new method works in both the pure absorption and pure scattering regimes, and properly treats anisotropic scattering processes. The method comprises a simple addition to the MC estimators required for radiative and photoionisation equilibrium calculations, and does not therefore represent a substantial computational overhead to the MC RHD methods described by Haworth & Harries (2012).

In this second paper of the series we apply the new method to the follow the formation of a massive star in three dimensions. In Section 2 we review our numerical method and in Section 3 we describe our initial conditions. In Section 4 we describe the evolution of the model, and in Section 5 we compare our simulation with observational constraints. We discuss our model in the context of recent observations in Section 6 and we make our conclusions and discuss future directions of our research in Section 7.

2 Numerical method

We employ the torus radiation hydrodynamics code (Harries, 2000; Harries et al., 2004; Harries, 2011, 2015) which uses a MC based treatment of the radiation field coupled with an adaptive mesh hydrodynamics scheme and a Lagrangian sink particle implementation. We direct the reader to the first paper in this series (Harries, 2015) for a detailed description of the numerical methods used to solve the RHD equations. We emphasise that we have extensively tested and benchmarked the MC RT scheme against dust continuum benchmarks (Harries et al., 2004; Pinte et al., 2009) and H ii region models (Haworth & Harries, 2012; Bisbas et al., 2015) as well as time-dependent heat-transfer benchmarks Harries (2011). We have also tested the coupled radiation-hydrodynamics scheme against the usual 1-D benchmarks, such as the Sod shock tube, the expansion of an H ii region (Harries, 2015; Haworth et al., 2015), and radiation-driven shell tests (Harries, 2015). The Lagrangian sink particle scheme has been benchmarked against standard n-body tests as well as ‘standard’ star formation tests such as Bondi accretion, Shu collapse, and Bondi-Hoyle accretion models (Harries, 2015).

3 Initial conditions

Our calculation starts with a 100  M⊙ cloud with a 0.1 pc radius and an $\rho(r)\propto r^{-2}$ density gradient. We impose solid body rotation with an angular frequency of $5\times 10^{-13}$  s*-1*. We assume a $10^{-12}$ g cc*-1* threshold density of accretion onto the sink particle (i.e. mass above this threshold within the sink accretion radius is assumed to accrete onto the sink). We start with a 0.001  M⊙ sink particle, which is held fixed at the origin. Our computational domain is an $8\times 10^{17}$ cm-sided cube, with a minimum of 7 and maximum of 12 levels of refinement. The grid is refined in nested spherical volumes, with the inner 150 au refined at 13 au resolution, then 26.1 au resolution beyond this out to 300 au, 52.2 au out to 500 au, 104 au out to 1000 au, 209 au to 2100 au, and 418 au for the rest of the domain. For the purposes of this initial calculation the grid is not refined or unrefined as the calculation proceeds. We impose outflow/no inflow conditions at all boundaries.

3.1 Dust model

We assume the dust is composed of silicate grains (Draine & Lee, 1984) in an MRN (Mathis et al., 1977) size distribution ($n(a)\propto a^{-3.5}$) with a minimum size of 0.005 $\mu$m and a maximum size of 0.25 $\mu$m. For the purposes of this calculation we keep this size distribution constant in time and space, although we note that detailed calculations suggest that coagulation and shattering may lead to significant modification of the dust size distribution during collapse (Suttner et al., 1999). We further assume that the gas and dust are strongly dynamically coupled, and that the dust sublimates when its temperature exceeds 1500 K.

3.2 The protostar model

We take the mass of the sink particle and interpolate in the $\log(\dot{M})=-3$ evolutionary models of Hosokawa & Omukai (2009) to obtain the luminosity and radius of the protostar. We use the sink particle accretion rate to calculate the accretion luminosity (using the Hosokawa & Omukai (2009) radius) and assume that all the accreted kinetic power is converted to protostellar luminosity. The spectral energy distribution is found by interpolating in a grid of Kurucz model atmospheres using the temperature and surface gravity of the protostar, and adding on a blackbody contribution for the accretion luminosity with a temperature derived by assuming the accretion luminosity is liberated uniformly across the protostellar surface.

4 Results

The accretion rate at early times is about $10^{-3}$ M⊙ yr*-1*, rising to a maximum of $2\times 10^{-3}$ M⊙ yr*-1after 6 kyr when the accretion rate starts to decline as the disc forms and effects of radiation feedback start to grow (the protostellar mass has reached about 10 M⊙ at this point, see Figure 1). Note that a non-rotating model with our mass and density profile that is ballistically infalling would have an accretion rate of $2\times 10^{-3}$  M⊙ yr-1* (Whitworth & Ward-Thompson, 2001). We plot the instantaneous accretion rate in Figure 1, which is the accretion rate at the time the code dumps the hydrodynamical data (every 100 yrs). The accretion rate show significant stochastic variability throughout the run, at early times due to the feedback between increases in accretion luminosity (which dominates at this stage) and growth in thermal pressure. The protostellar growth is characterised by stochastic disc accretion at later times.

In Figures 2 and 3 we show the density and velocity distributions of the model via slices through the $y$-axis and $z$-axis respectively. Radiation pressure starts to produce outflow cavities at around 10 kyr, but because the disc is not yet fully established the outflow on one side of the disc dominates. The system because more symmetrical with time as two radiation-pressure dominated, high velocity cavities emerge with a large, thick, accretion disc surrounding the protostar. The cuts through the $z$-axis show the emergence of a large (1500 au radius) accretion disc, within which are large scale spiral structures driven by gravitational instabilities.

5 Observational comparisons

How well does our model agree with observational constraints on massive star formation? We can make such comparisons at two levels. The first is to compare the gross characteristics (masses, accretion rates, luminosities) from our simulation with those found from observation. The second is to make a comparison at a level that is closer to the observations themselves, by synthesising how the model would look to an observer via so-called simulated observations.

5.1 Accretion rate

We start by comparing the accretion rate from the simulation with observations. In a recent review of discs around luminous young stellar objects Beltrán & de Wit (2016) presented compilations of observational data of accretion rates as a function of stellar mass (their Figure 15). They explored four different estimators of the accretion rate: (i) The infall rate estimated from red-shifted absorption line profiles, (ii) by assuming that the infall velocity is equal to the rotational velocity, (iii) by assuming that all the material from the rotating structure around the protostar will accrete on a free-fall time, and (iv) by assuming that the accretion rate is one sixth of the outflow rate derived from observations of molecular outflows.

In Figure 4 we show the evolutionary track of our model in the $\dot{M}-M_{*}$ plane against the data collated by Beltrán & de Wit (2016). We show two different measures of $\dot{M}$ from our simulation, the first being the accretion rate onto the protostar (smoothed over 100 yr intervals), and the second being the mass flux passing through a spherical surface of radius 1500 au (the approximate radius of the protostellar disc, see below). The latter rate might be regarded as more representative of the observational rates, which probe large scales. The agreement with observation is generally good (although there is significant scatter in the data), but the model accretion rates seem to be systematically at the lower end of the distribution of observational rates (apart from those derived by scaling the outflow rates).

At the end of our simulation the bolometric luminosity of the system is $1.77\times 10^{5}$ L⊙, of which the accretion luminosity contributes $3.7\times 10^{4}$ L⊙ (giving a ratio $L_{\rm acc}/L_{\rm bol}=0.2$). The values are in good agreement with the high mass sample presented in Beltrán & de Wit (2016).

5.2 Disc mass and radius

We work on the basis that the disc in the numerical simulations corresponds to the volume within which the gas mass density exceeds $10^{-15}$ g cc*-1*. Although apparently arbitrary, this threshold matches a by-eye estimate of the disc extent from visualisations of the simulations (see Figure 3), and we note that the same density threshold was employed by Krumholz et al. (2009) when examining the properties of their model disc. The disc is not completely circular at this density cut due to large spiral density waves, but varies between 1000 and 2000 au in radius, and we adopt 1500 au as the characteristic radius of the disc. We find the disc mass by summing the density of all cells within 1500 au of the protostar that has a density exceeding $10^{-15}$ g cc*-1* . By the end of the simulation the disc mass is about 7  M⊙, or about 30 percent of the protostellar mass.

We compare our disc mass and radius with observationally derived qualities in Table 1. Broadly speaking the gross properties of our model disc are comparable with observations, falling into approximately the middle of the distribution in terms of mass, radius, and disc-to-protostellar mass ratio.

We extracted the $y$-component of velocity of cells in our model disc which have densities greater than $10^{-15}$ g cc*-1* and that have a face in the $z=0$ plane (i.e. the mid plane of the disc). We find that the disc is in Keplerian rotation (see Figure 5). We note that Ilee et al. (2016), Johnston et al. (2015) and Chen et al. (2016) identify Keplerian rotation in their observations of massive discs. However the observational measurements of the disc kinematics are complicated by the fact that they are obtained from molecular line tracers that may also simultaneously probe the rotating and infalling envelopes (for example Hunter et al. (2014) observed a combination of Keplerian rotation and infall). In order to make a more realistic comparison between the model and observations it is necessary to compute synthetic observations from the model and test the model more directly via molecular line profiles (see Section 5.5).

5.3 Spectral energy distributions

We calculated the spectral energy distributions (SEDs) from the models density and temperature distributions assuming an inclination of 60∘, and these are plotted in Figure 6. Strong 10 $\mu$m absorption features are seen in all spectra, consistent with its deeply embedded nature and the presence of small silicate grains. The rapid growth of the near- to mid-IR flux is seen as the model evolves and the central object becomes more luminous.

Since we know that the model is meant to represent a MYSO it is worth testing whether our object would appear in the Red MSX Survey catalogue (Lumsden et al., 2013), which collates data on a large number of MYSOs selected from the point source catalogue of the MSX satellite (Price et al., 2001). In order to do this we constructed synthetic fluxes in the appropriate wavebands assuming a canonical distance of 1 kpc. The cuts used by Lumsden et al. (2013) are $F_{21}>2F_{8}$, $F_{21}>F_{14}$, $F_{14}>F_{8}$, $F_{8}>5F_{K}$ and $F_{K}>2F_{J}$. (Note that since our model is still deeply embedded at the end of the simulation we neglect the last constraint.). We plot the colour cuts in Figure 7, and do indeed find that the first three cuts are easily passed at all times. However we note that the $F_{8}>5F_{K}$ constraint is not passed for the later models, and at we attribute this to a deficit of $F_{8}$ flux due to the presence of a strong silicate 10 $\mu$m feature.

5.4 Images

Before we proceed to examine synthetic images of our models it is worth taking a short digression into the jargon of the field of massive star formation, focussing on the term ‘bipolar outflow cavities’. At early times our simulations show (see Figure 2), a pair of roughly diametrically opposed rarified regions that are formed by a rapid ($\sim 100$  km s*-1*) quite collimated, outflows driven by radiation pressure. As the simulation progresses and the disc forms, the outflows straighten towards the rotational axis and the hot gas inside the rarified regions expands to form a pair of ‘bubbles’. However, bipolar cavities may also be carved by highly-collimated magnetically-driven outflows (jets), and these are typically characterised by much smaller opening angles. The self-consistent calculation of the formation of these magnetic towers requires full radiation-magneto-hydrodyamics at a resolution that would render large-scale simulations intractable, but some assessment of their role in modifying the effects of radiation pressure may be made by incorporating a sub-grid model of the bipolar outflow e.g. Cunningham et al. (2011) and Federrath et al. (2014). In the following discussion we refer to bipolar cavities, and cavity walls, in reference to the rarified structures that are formed along the rotation axis of our simulations, and defer an investigation of magnetically-driven outflows to future work.

We calculated synthetic images of the model at 35 kyr. These monochromatic images were calculated at the nominal wavelengths of the MSX passbands and also Herschel PACS 70 µm, Herschel SPIRE 250 µm and 350 µm, and ALMA 850 µm. The 2 µm image is dominated by protostellar light scattered in the walls of the bipolar cavities. The mid-IR images show significant thermal emission from the warm dust at the base of cavities, and the 70 µm is produced by dust throughout the cavities and the cavity walls. The longer wavelength images are dominated by cool, dense dust in the disc around the protostar.

Alvarez et al. (2004) presented H and K band speckle imaging of outflows in MYSOs. They found that 6 objects from their sample of 21 showed a conical nebula which could be interpreted as arising from scattered light from the dusty walls of the outflow cavity. Typically these nebulae had an extent of about half an arc second, and are monopolar in nature, the implication being that scattering from the more distant cavity is obscured by intervening material such as a disc. Some of the objects were found to have a gap between the location of the protostar and the base of the nebular cone. Our 2 µm image (top left of Figure 8) is consistent with the Alvarez et al. (2004) observations, although we note at 60∘ inclination the image peak occurs at the position of the protostar since there is substantial unscattered near-IR thermal emission from the unresolved inner disc and protostar. However, although there is some low surface brightness emission at large angular distances from the protostar, the image is dominated by a conical, monopolar nebula whose centre of brightness is offset from the location of the protostar. The angular extent of the conical nebula (about 2 arc seconds) is rather larger than the typical size of those observed by Alvarez et al. (2004) at our 1 kpc distance.

de Wit et al. (2009) obtained diffraction limited (0.6 arc second) 24.5 µm images of 14 massive star formation regions. These images revealed the presence of MYSO sources that are resolved on arc second scales, and that they are (to first order) circular on the sky. They found that it was possible to simultaneously fit the imaging and SED data for the objects using spherical models with a $p=1$ density profile, although there were a subset of objects which were not well-fitted by spherical models. They conclude that these objects are viewed either edge-on or face-on. Our 21 µm image (top right of Figure 8) is dominated by thermal emission from the disc and the warm cavity walls. We note that the peak surface brightness of our model (about 107 MJy/str) is similar to those found by de Wit et al. (2009) but clearly the image is not circularly symmetric on the sky.

Olguin et al. (2015) found that three MYSOs at distances between 4–5 pc were marginally resolved at 70 µm in Herschel Hi-GAL survey data, and also found that the data could not be adequately fit with spherical models. Using two-dimensional RT models it was found that a lower power-law density index (0.5) than that of pure infall (1.5) was necessary to fit the data, and they attributed this to extended emission from warm cavity walls. We re-scaled our fiducial 70 µm image to 5 kpc distance and smoothed it to Herschel PACS resolution. We find that our model image would be very marginally resolved at this distance and resolution, but that an observer might be able to discern that the object was elongated.

As mentioned in Section 4 the disc shows significant spiral structure. In order to ascertain whether this predicted structure could be observed by the Atacama Large Millimetre Array (ALMA) we constructed a pole-on monochromatic thermal dust continuum image at a wavelength of 1 mm using the 35 kyr model. This was then used as the basis of a synthetic ALMA observation using CASA (McMullin et al., 2007). For this image baselines were combined from extended and intermediate ALMA configurations (finished ALMA configurations 6 and 20, as included in the CASA software package) and the total time used was 80 minutes (1 hour for the extended configuration and 20 minutes for the intermediate). Thermal noise was added to the image using an assumed amount of precipitable water vapour of 0.5 mm. The resulting image was cleaned interactively down to the level of 0.2 mJy/beam. The resulting synthetic ALMA images is shown in Figure 9 (right panel) with the radiative-transfer simulation used as the skymodel for comparison (left panel). One can see that the large-scale structure observed in the radiative-transfer model is beautifully recovered in the CASA simulated ALMA observation. Furthermore the sublimation radius (seen as the black spot in the centre of the radiative-transfer model) is also observed in the CASA simulation as a dip in flux in the centre of the model.

5.5 Molecular tracers

Many of the recent advances on the dynamics of discs around MYSOs have been dependent on high resolution molecular observations. We have not treated time-dependent chemistry during our simulation, but by assuming a simple prescription for the abundance distribution of the target molecule we can use the molecular line module of torus to simulate the line emission data cube based on the RHD model’s density and temperature distribution.

We focus on transitions of the methyl cyanide molecule (acetonitrile, CH3CN) since the high excitations transitions of the $k$-ladder have been used to probe massive star discs. We adopt the same abundance distribution relative to H2 as Johnston et al. (2015) i.e.

[TABLE]

Using this prescription we found that the emission from the system had a significant component from the hot walls of the bipolar cavities and very strong emission near the sink particle, which is not something that is observed. We therefore added the further constraint that only gas with a mass density of greater than $10^{-15}$ g cc*-1* and a gas temperature of less than 1500 K had a non-zero CH3CN abundance.

The molecular line calculations were computed with torus, using the molecular RT physics module as described and tested in Rundle et al. (2010). We used the CH3CN energy levels from the lambda database (Schöier et al., 2005) and we took the collisional rates from Green (1986). We focus on the $J=13\rightarrow 12$, $K=3$ (239.096411 GHz) transition, which has a critical density of $4\times 10^{6}$ cm*-3* at 500 K. The gas density in the protostellar disc exceeds $10^{8}$ cm*-3* and we therefore assumed local thermodynamic equilibrium, adopting the gas temperatures from the RHD calculation. Using the $t=35$ kyr RHD snapshot we computed a data cube with $256\times 256$ spatial pixels (corresponding to $6.2$ arc seconds on a side at an assumed distance of 1 kpc) and 33 velocity channels ranging from $-8$  km s*-1* to $+8$  km s*-1*. An inclination of 30∘ was adopted.

Figure 10 shows the total integrated spectrum of the $J=13\rightarrow 12$, $K=3$ transition (upper left). A double-peaked profile is seen, with a significantly stronger blue peak due to the asymmetric nature of the disc. The brighter blue peak is highly dependent upon the viewing angle to the model however, a different choice of viewing angle can result in a stronger red peak or two equally bright peaks. The continuum-subracted zeroth order moment map (upper right) shows that the emission is strongly peaked towards the centre of the disc, as expected for this high excitation line. The integrated intensity map also shows a gap in emission at the very centre of the disc where CH3CN is depleted due to the high temperature. The first-order moment map (lower left) shows the expected pattern for a Keplerian disc. The position velocity diagram (lower right) shows that the rotation curve of the disc could be extracted from these CH3CN observations if they were made with sufficient resolution and sensitivity. We note that similar calculations were performed by Krumholz et al. (2007), and although they assumed a constant CH3CN abundance, they drew broadly the same conclusions regarding the feasibility of extracting kinematic information from the observations.

Four of the five observations listed in Table 1 (the observation of IRAS 16547 does not include CH3CN as one of its targets) included measurements of a methyl cyanide (CH3CN) line, and in each of them the velocity extent is 10–15 km s*-1* across the line. This is in agreement with the $\sim$10 km s*-1* seen in our synthetic observations of the CH3CN $K=3,J=13\rightarrow 12$ line. The peak surface brightness in this transition is slightly stronger (a factor of 1.5–2) than the observations of Chen et al. (2016) and Ilee et al. (2016).

6 Discussion

Our simulation is consistent with a picture of massive star formation as a scaled-up version of the low-mass paradigm, with stochastic disc accretion onto a protostar from a Keplerian disc accompanied by bipolar outflows. Since we have adopted a novel radiative-transfer scheme for this calculation it is interesting to compare our results to other recent studies. Klassen et al. (2016) studied high mass star formation by implementing the hybrid radiation transport scheme in the flash AMR hydrodynamics code. They started with spherical protostellar cores of 30 M⊙, 100 M⊙, and 200 M⊙, with an $\rho(r)\propto r^{-3/2}$ density profile and the same solid body rotation as this study. Broadly the result of their 100 M⊙ simulation is very similar to our calculation, with the presence of a thick disc showing evidence of spiral structures and fast, bipolar cavities that initially show significant asymmetries (these authors also attribute these asymmetries to resolution effects due to the carteisan mesh). The accretion rate onto the star is similar to ours during the early phases of the calculation, and both show the accretion rate initially peaking after 5–10 kyr. However the Klassen et al. (2016) shows a rapid growth in the accretion rate as the disc becomes gravitationally unstable. We attribute the difference to the fact we have fixed our protostar at the origin for the purpose of this first calculation, which will suppress gravitational instabilities that would otherwise enhance accretion. Despite this difference the similarity between the density slices in Figure 2 and those of Figure 13 in Klassen et al. (2016) is remarkable.

There has been considerable controversy in the literature regarding the role of radiatively-driven Rayleigh-Taylor (RRT) instabilities during the formation of massive stars. First suggested by Krumholz et al. (2009) RRT instabilities arise from when the dense gas in the cavity shells starts to penetrate the rarified gas of the cavities via finger-like intrusions. Numerical simulations focussing on galactic scales confirm the presence of such instabilities in optically-thick radiation-pressure driven dusty winds (Krumholz & Thompson, 2012; Davis et al., 2014; Rosdahl & Teyssier, 2015; Tsang & Milosavljević, 2015), although the details of the dynamics are found to be sensitive to the adopted radiation treatment.

Kuiper et al. (2012) argued that the RRT instabilities were a result of the underestimate of the radiative driving due to the use of grey flux-limited diffusion by Krumholz et al. (2009), whereas the hybrid method they employed better captured the momentum deposited by the radiation at the $\tau=1$ surface. It appeared that the results of Klassen et al. (2016), who also used a hybrid method, supported this finding, since they did not observe RT instabilities during their simulations. Recently Rosen et al. (2016) returned to this issue and implemented a hybrid method into the same hydrodynamics code used by Krumholz et al. (2009). They found that although the implementation of the hybrid scheme suppresses the onset of the RT instabilities they nonetheless still occur. They argue that since the star is fixed in the Kuiper et al. simulations the optical depth effects that can stochastically occult the radiation field and suppress the radiation driving are absent, leading to a more stable outflow. They also argue that the Kuiper calculations (which are performed on a fixed spherical mesh) do not have a sufficient resolution at the cavity boundary ($\sim 400$ au) to capture the growth of the small scale perturbations that should seed the instability. They point out that although the Klassen simulations use an adaptive mesh they do not refine on the gradient of the radiation density, and thus the resolution along the cavity ($\sim 160$ au) is insufficient to resolve the instability. Note that for comparison the calculation presented here has a resolution of 418 au at a radius from the protostar of 4000 au (roughly the radius at which the instabilities occur in the Rosen et al. simulation).

The question seems to be not whether RRT instabilities can occur, but whether a particular simulation is conducted in a regime where the instability is important. It seems likely that in some published simulations the radiation pressure of the direct (protostellar) radiation field at the dust sublimation radius imparts sufficient momentum that the outflow cannot be halted by infalling material in the outer envelope, and reprocessed radiation plays no role in the driving process. Such a simulation would not be expected to show RRT instabilities. However for different initial conditions the momentum imparted to the gas at dust sublimation radius may be insufficient to drive a fast, stable outflow and instead the reprocessed radiation field is an important driving mechanism, and RRT instabilities result. Of course the preliminary simulation presented here neither has neither a moving star nor a grid that refines on the radiation field, so although we find good agreement with the Klassen et al. (2016) simulations we defer a detailed discussion of the role of RT instabilities to future work.

Our RHD model finishes with a luminous protostar, surrounded by a thick disc, large cavities, and an optically think infalling envelope. The bolometric luminosity, and the fraction of accretion luminosity to protostellar luminosity are in good agreement with observations. The model SEDS (see Figure 6) are comparable to those in the RMS survey (see, for example, Mottram et al. 2011). Of course for these deeply embedded objects much of the leverage of the SED on the density structure is lost as the radiation field is reprocessed by the optically thick envelope prior to being observed.

Extending the comparison to imaging we predict near-IR monopolar conical nebula caused protostellar and inner-disc emission scattering off the cavity walls. The extent of the scattered light nebula appears rather larger than the observations, although the morphology of the nebula is sensitive to the adopted viewing angle. The presence of the bipolar cavities also has a strong impact on the mid-IR images, with the 21 µm image particularly showing extended emission from the cavity walls. The fact that 24.5 µm imaging shows mainly circular symmetric emission places a strong constraint on our models, and suggests that our outflow cavities are overly large in the RHD model.

It is interesting to compare our synthetic images with those from analytical prescriptions of the disc and envelope density structure published in a sequence of papers by Zhang and co-workers (Zhang & Tan, 2011; Zhang et al., 2013, 2014). The outflow cavity in these models is assumed to arise from a disc wind which evolves in such a way that the cavity opening angle increases with time, reaching approximately 30∘ when the protostar is 10  M⊙ and extending to 90∘ when the protostar is 24  M⊙. Unfortunately the analytical density distribution for the 24  M⊙ model in Zhang et al. (2014; hereinafter Z14) is sufficiently different to our RHD simulation that the images do not compare well. However, the 16  M⊙ model images (Figure 21 of Z14) do show some obvious qualitative similarities with our synthetic images (Figure 8). In particular the 20 µm image show is significantly extended in Z14, in agreement with our models and in contradiction with observations. The 2 µm scattered-light image in Z14 is even more extended than our own, and is therefore also large compared to observations.

The root cause of this disagreement is probably the lack of collimated magnetically-driven jets in our calculations (and those of Z14), which may result in bipolar cavities with smaller opening angles at later times. We are investigating this using two-dimensional simulations and the sub-grid jet model of Federrath et al. (2014), and will report on these simulations in a future paper.

Fortunately as we go to longer wavelengths the agreement with observation improves, and we predict that our system would only be marginally resolved at 70 µm by Herschel PACS, in agreement with the findings of Olguin et al. (2015). Processing a 1 mm dust continuum image with the CASA software suggests that ALMA should be able to resolve the spiral waves that occur in our disc, and even perhaps the dust sublimation radius.

Our synthetic molecular line observations suggest that kinematical data should be straightforwardly extracted (see Figure 10), allowing a straightforward determination of the protostellar mass by assuming the disc is in Keplerian rotation. This is a good assumption for our RHD model (see Figure 5), but the observations paint a more complex picture, in which the infalling envelope has an effect on the molecular line RT. Clearly the main limitation is the lack of a chemistry model to predict the abundance distribution of CH3CN. Although it is possible to couple RHD with time-dependent chemistry for limited reaction networks it would currently be intractable to attempt such a computation of a network detailed enough to follow the formation of CH3CN. However an equilibrium chemistry post-processing of the RHD model is currently within our reach, and we are working on implementing this with the torus framework.

7 Conclusions

We have presented a new simulation of the formation of a massive star using a novel MC-based RHD method. Our model shares characteristics with RHD simulations based on hybrid RT methods (e.g. Kuiper & Klessen 2013, Klassen et al. 2016), such as a large, thick protostellar disc and fast, broad, outflow cavities. The bulk of the protostellar mass is acquired via stochastic disc accretion.

We do not observe radiatively-driven Rayleigh-Taylor instabilties, but as noted above we do not refine our AMR mesh at the cavity shell, which may mean that we may under-resolve the small-scale variations that seed the instability (Rosen et al., 2016). However we are currently in the process of running a suite of 2D (cylindrical) models and we find that at early times the bipolar outflows may be overwhelmed by infalling material. This is primarily due to variability in the protostellar accretion rate changing the driving luminosity, but it sensitive to initial conditions. In order to fully resolve this issue we are implementing a grid refinement method that refines the AMR method at the cavity boundary.

We have attempted to compare our model with observations. We find that our model is consistent with observations in terms of protostellar luminosity and accretion rates, and that the protostellar disc matches the observed distribution in terms of disc mass and radius. We predict the disc to be in Keplerian rotation, but with significant large-scale gravitationally induced spiral structure that would resolvable with ALMA.

The SEDs predicted by our models show good agreement with observed MYSO spectra, and that the synthetic colours of the object would mean that it would be selected as an MYSO using most of the colour cuts adopted for the RMS survey. However, we note that the strength of the 10 µm silicate feature may be sufficient to suppress the 8 µm flux in such a way that it might fail the $F_{8}>5F_{K}$ cut.

Synthetic images of RHD model show some features that are observed, such as conical scattered light nebula in the near-IR and 70 µm emission that would only be marginally resolved in Herschel PACS data. Although the peak surface brightness of our synthetic 21 µm imaging matches the observations well, its angular extent does not, and this perhaps indicates that our RHD model is over-predicting the extent of the bipolar cavities.

We show that synthetic molecular line observations of CH3CN predicts the correct surface brightness in the line, as well as the line velocity width. We find a blue-red asymmetry in the line profile, and attribute this to asymmetries in the disc itself (as opposed to RT effects as the line emission from the disc passes through the infalling envelope). We note that our molecular line calculations are limited by a lack of self-consistent chemistry.

In summary, our RHD model shows similar characteristics to other contemporary studies in terms of the disc, envelope and outflow properties. We have shown that these properties are broadly consistent with observations across a wide variety of observables. In future studies will address the issues of radiatively-driven Rayleigh Taylor instabilities, and study the disc fragmentation and binary formation. We will address the effect of including a sub-grid model of the collimated magnetically-driven outflow on the accretion process. We also hope to couple the torus code with a chemical model in order to make more self-consistent predictions of molecular line diagnostics.

Acknowledgments

The calculations for this paper were performed on the University of Exeter Supercomputer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital fund of BIS, and the University of Exeter, and on the Complexity DiRAC Facility jointly funded by STFC and the Large Facilities Capital Fund of BIS. TJH and TAD acknowledge funding from Exeter’s STFC Consolidated Grant (ST/M00127X/1). We thank Takashi Hosokawa for kindly providing us with the protostellar evolutionary model. We are grateful to Maite Beltran for providing the data for Figure 4, and we thank Dave Acreman, John Ilee and Tom Haworth for useful discussions. We thank the anonymous referee for a helpful report.