© ESO 2020

1 Introduction

During the formation of massive stars (≳8 M_⊙), radiation pressure becomes important against the gravity of the collapsing molecular cloud. The formation of an accretion disk with polar outflows provides a mechanism for circumventing the radiation pressure barrier (see, e.g., Kuiper et al. 2010a) and allows the forming star to become massive. This disk is expected to fragment and produce companion stars. There is growing observational evidence that supports this scenario. Observations of disks around massive (proto-)stars are reported by, for example Johnston et al. (2015), Ilee et al. (2016), Cesaroni et al. (2017), Ginsburg et al. (2018), and Maud et al. (2019). Some of these disks have also been shown to be Keplerian-like.

Moreover, there is evidence that, early in evolution, these disks gain enough mass to become self-gravitating, form spiral arms, and fragment. Ilee et al. (2018) observed a fragmented Keplerian disk around the proto-O star G11.92-0.61 MM1a, with a fragment MM1b in the outskirts of the disk, at ~2000 au from the primary. Beuther et al. (2017) reported a smaller disk-like structure around the central object in the G351.77-0.54 high-mass hot core, and a fragment at about ~1000 au. Johnston et al. (2020) also observed spiral arms and instability in a disk of radius ~1000 au around the O-type star AFGL 4176 mm1. Maud et al. (2019) reported substructures in a Keplerian disk around the O-type (proto-)star G17.64+0.16.

Several studies have used three-dimensional hydrodynamical simulations of a collapsing cloud, including radiation transport and self gravity, in order to study the formation of a massive star. Some of them lead to a massive star surrounded by a fragmenting accretion disk, and some form stars via filament accretion (we offer a literature review, including methods, initial conditions, and results in Sect. 9). Krumholz et al. (2009), Kratter et al. (2010), Rosen et al. (2016), Klassen et al. (2016), Girichidis et al. (2012), and others used adaptive mesh refinement (AMR) grids and sub-grid sink particle models. In the case of Girichidis et al. (2012), as many as approximately 400 sink particles were reported in some of their runs, or in the case of Rosen et al. (2016), up to about 30; in contrast, Klassen et al. (2016) report none. These large disparities in the number of possible formed companions under similar conditions provokes questions on the role played by spatial resolution and the sink particle algorithms used in these studies in the final outcome of the system. In Sects. 8 and 9, we explore this matter in detail, and find that higher resolution is needed to resolve the Jeans length than used in previous studies.

Accretion in fragmented disks is expected not to be a smooth process, but to be characterized by episodic accretion of some fragments. The release of gravitational potential energy creates an increase in luminosity, i.e., an accretion burst. Accretion bursts offer an explanation for the luminosity burst events in regions of massive star formation which have been reported by Hunter et al. (2017), Caratti o Garatti et al. (2017), and Sugiyama et al. (2019) for example. Recently, Chen et al. (2020) reported the observation of disk substructures associated with an accretion burst event, thus providing a link between thetwo phenomena. The simulations in Meyer et al. (2017, 2018, 2019) show accretion bursts; these latter authors performed an extensive analysis of the accretion process and the intensity and frequency of the bursts.

Massive stars do not typically form in isolation. Studies on multiplicity in more evolved systems of massive stars show that a large portion of them are spectroscopic binaries (see, e.g., Sana et al. 2012; Kobulnicky et al. 2014; Dunstall et al. 2015). The observationsof fragmenting accretion disks suggest the possibility that close binaries may form by disk fragmentation and later inward migration of the companions (Meyer et al. 2018).

In this study, we present the highest-resolution self-gravity-radiation-hydrodynamical simulations of disk fragmentation around a forming massive (proto)star to date. We continue the approach taken by Meyer et al. (2018), in which no sink particles are set purposely, so that fragmentation is described self-consistently as the interplay of self-gravity of the gas, cooling of the gas, shear of the disk, and other gravito-radiation-hydrodynamical processes. Instead, we used a time-independent grid whose spatial resolution scales logarithmically with radius, allowing very high resolution in the areas of most interest. We set up a series of simulations with increasing resolutions, and explore the physical processes that occur in fragmentation (how fragments evolve and interact), the properties of the fragments (orbits, mass, and temperatures), their fate (whether they evolve further and eventually become companion stars), and the effects that resolution has on fragment statistics.

The outcome of our simulation was already used in Ahmadi et al. (2019) to produce synthetic observations with ALMA and NOEMA. In that paper, many observationally relevant quantities were computed, including 1.37 mm continuum images, integrated intensity and intensity-weighted peak velocity maps of CH₃ CN(12₄-11₄), position–velocity plots, and Toomre Q maps, taking into account several inclinations. We refer the reader to that article for an ample discussion on the observability of the results presented here.

This paper is organized as follows. Sections 2 and 4 present the methodology of both the simulations, and the post-processing of the generated data. Sections 3 and 5–7 present our results and a discussion of disk fragmentation, the fragments themselves, and companion formation. Section 8 contains a convergence study of our results, and Sect. 9 a literature review and comparison of our results with previous studies.

2 Setup of the simulations

2.1 Physics

In order to model the problem, we consider the three-dimensional collapse of an initially axially symmetric molecular cloud of total mass 200 M_⊙ and radius 0.1 pc. The initial density and rotational profiles are given in Sect. 2.3.

The cloud dynamics are modeled as an ideal gas with the hydrodynamics equations: $$\begin{array}{l}\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho v)=0,\\ \frac{\partial }{\partial t}(\rho v)+\nabla \cdot (\rho v\otimes v+PI)=\rho a_{\text{ext}},\\ \frac{\partial E}{\partial t}+\nabla \cdot \left((E+P)v\right)=\rho v\cdot a_{\text{ext}},\end{array}$$

where ρ, P, and E are the density, pressure, and energy density, respectively; v is the velocity, and a_ext is the acceleration source term. These equations are solved with the hydrodynamics module of the numerical grid code Pluto (Mignone et al. 2007), with additional modules for handling radiation transport and self-gravity.

The gravity of the forming massive star, the self-gravity of the fluid, and the radiation forces are incorporated into the code via the acceleration source term,such that a_ext =a_⋆ +a_sg +a_rad, and $$\begin{array}{l}{a}_{\star }=-\frac{G{M}_{\star }}{r^2}{e}_r,\\ a_{\text{sg}}=-\nabla {\Phi }_{\text{sg}},\text{where}{\nabla }^2{\Phi }_{\text{sg}}=4\pi G\rho ,\\ a_{\text{rad}}=-\frac{\nabla \cdot F_{\star }}{\rho c}{e}_r-{\kappa }_R\frac{D\nabla E_R}{c}.\end{array}$$

The variables involved in a_rad are detailed below. The Poisson equation of self gravity is solved by means of a diffusion ansatz. More details on the implementation of the Poisson solver for the self-gravity module are given in Kuiper et al. (2010a).

The gas is assumed to be calorically perfect, that is, governed by the calorical equation of state P = (γ − 1)E_int (where γ = 5∕3, the value for H₂ at low temperatures, and E_int is the internal energy density).

Radiation transport is incorporated by the method described in Kuiper et al. (2010b), but advanced to a two-temperature approach, which we summarize in the remainder of this section. We treat the star as an emitter of radiation, which the dust and gas present in the cloud absorb and re-emit in addition to the thermal radiation due purely to hydrodynamic compression. This treatment allows us to divide the total radiation flux F_tot into two contributions: the flux from the star, F_⋆, and the flux from thermal (re-)emission from dust and gas, F.

Stellar irradiation is solved by means of the time-independent radiation transport equation (with the valid assumption that the photon travel time is small compared to the time-step of the simulation, and only considering absorption), which yields $$F_{\star }(\nu ,r)=F_{\star }(\nu ,R_{\star }){\left(\frac{R_{\star }}{r}\right)}^2{e}^{-\tau (\nu ,r)},\text{where}\tau (r,\nu )={\displaystyle {\int }_{R_{\star }}^r\kappa }(\nu ,r)\rho \text{d}r$$(7)

is the optical depth, κ(ν, r) the frequency-dependent opacity, and R_⋆ the stellar radius. This equation is integrated for each ray direction and frequency bin in each time-step of the hydrodynamical simulation.

The frequency-dependent opacity κ(ν) is the sum of the dust and gas opacities. For the dust opacities, we use the tabulated values by Laor & Draine (1993), including 79 frequency bins, and we calculate the local evaporation temperature of the dust grains using the formula of Isella & Natta (2005). The gas opacity is set to a constant value of 0.01 cm² g⁻¹. We also require the flux of the forming star at the stellar radius R_⋆, for which we use tabulated evolutionary tracks for accreting high-mass stars calculated by Hosokawa & Omukai (2009).

For the remaining thermal (re-)emission contribution, F is assumed to be frequency-independent. First, we take the zeroth moment of the radiation transport equation $$\frac{\partial E_{\text{rad}}}{\partial t}-\nabla \cdot F=\rho {\kappa }_P(4\pi B_{\text{rad}}-c{E}_{\text{rad}}),$$(8)

where B_rad = aT⁴ is the black-body energy density, c is the speed of light in vacuum, and κ_P is the Planck mean opacity. This equation is solved by means of the flux-limited diffusion (FLD) approximation, which consists in setting F = −D∇E_rad, that is, assuming that radiation transport can be treated as a diffusion problem. The diffusion constant is D = λc∕ρκ_R, where λ is the flux limiter, and κ_R is the Rosseland mean opacity. The change of internal energy of the system is $$\frac{\partial E_{\text{int}}}{\partial t}=-\rho {\kappa }_{\text{P}}(4\pi B_{\text{rad}}-c{E}_{\text{rad}})-\nabla \cdot F_{\star },$$(9)

where, for the ideal gas, E_int = c_vρT. We solve Eqs. (8) and (9) for the unknowns E_rad and T using the so-called two-temperature linearization approach, described in Commerçon et al. (2011a).

2.2 Geometry

The equations described in the previous section are solved in a three-dimensional, spherical, time-independent grid at five different resolutions, which are detailed in Table 1. The computational domain extends from r = 30 au to r = 20 626.5 au (= 0.1 pc). The coordinate grid is built as follows: the radial coordinate increases logarithmically from the inner boundary; the polar angle varies with the cosine function so that maximum resolution is achieved in the midplane (z = 0); and the azimuth is uniformly discretized. In the midplane the cells are approximately cubical, with two example cell sizes given as Δx in Table 1.

This choice of coordinates is explained and justified in more detail in Kuiper et al. (2011); briefly, a spherical grid guarantees strict angular momentum conservation around the central massive protostar, and the logarithmic grid in r allows for a focus on the phenomena closer to the forming massive star, while saving computational power. An explanation of how well this choice of grids resolves the phenomena studied here is presented in Sect. 9.2 together with a comparison to the grid choices of previous studies.

2.3 Initial and boundary conditions

The initial density profile is spherically symmetric, and has the general form $$\rho (r,t=0)={\rho }_0{\left(\frac{r}{r_0}\right)}^{{\beta }_{\rho }}.$$(10)

We choose β_ρ = −3∕2, and from the total mass and radius of the cloud, ρ₀ is determined to be ≈4.79 × 10⁻¹² g cm⁻³ at r₀ = 1 au. The value of β_ρ was chosen based on the results of previous simulations and observations of massive dense cores (see the discussion in Sect. 9.3).

The initial angular velocity is given by the profile $$\Omega (R,t=0)={\Omega }_0{\left(\frac{R}{R_0}\right)}^{{\beta }_{\Omega }},$$(11)

where R is the cylindrical radius, and we choose β_Ω = −3∕4. The ratio of kinetic energy to gravitational potential energy is set to 5%, which fixes the normalization parameter Ω₀ to ≈ 9.84 × 10⁻¹¹ s⁻¹ at R₀ =10 au. The normalization process is described in more detail in Meyer et al. (2018), together with a parameter scan for other choices of β_Ω. The selection of β_Ω = β_ρ∕2 keeps the ratio of kinetic energy to gravitational energy independent of the radius of the cloud; the value for this ratio in turn was chosen in accordance with the typical values found in Goodman et al. (1993) and Palau et al. (2014), for example. The initial radial and polar velocities are set to zero.

In summary, we give a low initial angular momentum to the cloud so that it forms a disk while collapsing. In a more realistic situation, inhomogeneities in the density profile may be present, as well as some initial subsonic turbulence and magnetic fields. However, the aim of this work is to study fragmentation of the disk due to the development of Toomre instabilities in the accretion disk that forms in the process, hence our choice of initial conditions. A more detailed discussion on the effects of the initial conditions and the inclusion of additional physical effects on the outcome of the simulations can be found in Sect. 9.3.

Both the inner and outer boundaries of the computational domain are semi-permeable, that is, matter is allowed to leave the computational domain but it cannot re-enter it. For the purposes of the calculations involving the forming massive star, all the mass in the sink cell is considered as accreted by the central massive protostar that forms within it. The initial temperature is 10 K throughout the computational domain, and the initial dust-to-gas mass ratio (used in the calculation of κ(ν)) is chosen to be 1% (Draine et al. 2007).

3 Overview of the temporal evolution of the system

In the following sections, we use the data from runs x8 and x16 as a reference for the description and analysis of our results. A discussion of the results for the other runs can be found in Sect. 8.

Figure 1 shows an overview of the time evolution of the system (for run x8) with maps of the density in the midplane (z = 0), the density in the xz-plane – to give an idea of the vertical structure of the system – the temperature, and the Toomre parameter (which is introduced in Sect. 4.1).

The cloud begins to collapse as soon as the simulation starts. Matter free-falls into the center of the cloud, where the sink cell is located, and eventually forms a massive star there. Simultaneously, angular momentum conservation yields an initially axially symmetric accretion disk that grows in size over time.

At around 4 kyr, when the disk is about 200 au in size, the axial symmetry is broken: spiral arms appear and then the disk fragments. The fragments form, interact with the disk, spiral arms, and other fragments, and can be destroyed. Very roughly, fragments have densities of more than ~ 10⁻¹¹ g cm⁻³ and temperatures higher than ~600 K, while spiral arms and other filamentary structures have densities higher than ~ 10⁻¹³ g cm⁻³. The background accretion disk has densities higher than 10⁻¹⁵ g cm⁻³. As shown inthe density maps of Fig. 1, the whole disk (including spiral arms) is relatively thin; the pressure scale height is shown below in Fig. 23. The spiral arms and filaments are dynamic: they continuously change shape, form, and merge in the course of an orbit, and they extend throughout the disk. Fragments are typically connected to each other and to the central massive protostar via spiral arms and filaments. Some fragments have the potential to form companion stars, as we describe in Sect. 7. Fragments can also be surrounded by secondary disks with their own spiral arms that can lead to the formation of more fragments. These mechanisms are described in Sect. 6.1.

The disk stabilizes at around 15 kyr due to the increase in mass of the central massive protostar and the effects of stellar irradiation in the innermost part of the disk, stopping fragmentation of the primary disk until the end of the simulated time. We refer to this period as the quiescent epoch. Some fragments survive the fragmentation epoch, but the majority of fragments migrate toward the central massive protostar or get accreted by it. Run x8 ends at t = 0.52t_ff, where $t_{\text{ff}}~\sqrt{R^3/(GM)}$ is the free fall time of the cloud; run x16 ends at t = 0.40t_ff.

Interestingly, the final state of run x8 is reminiscent of the young protostellar object G11.92–0.61 MM 1 reported in Ilee et al. (2018). This similarity was not intended a priori in the simulation setup. We compare the outcome of run x8 and the properties of the fragment produced to this potential observational counterpart in Sect. 6.8.

The central massive protostar formed by accretion from the disk is about 12 M_⊙ after ~ 14 kyr of evolution, as shown in Fig. 2a. The total luminosity of the protostar is calculated by summing the luminosity predicted by the evolutionary tracks of Hosokawa & Omukai (2009), and the accretion luminosity, that is, the gravitational energy released by accretion, and computed with L_acc = GM_⋆Ṁ∕(2R_⋆), where Ṁ is the accretion rate onto the (proto)star (Fig. 2c). The total luminosity is shown in Fig. 2b together with the stellar luminosity component. At around 12 kyr, the massive protostar starts burning hydrogen, and the stellar luminosity becomes comparable and later dominates over the accretion luminosity. Prior to that, the luminosity of the central massive protostar is dominated by accretion, which causes accretion-driven bursts; these can be seen in Fig. 2b and are discussed in Sects. 6.6.2 and 7.3.1.

Fig. 1 Time evolution of the system (run x8). In the elaboration of the plots, snapshots of the simulation were taken specifically at 3.0, 8.0, and 17.5 kyr for each respective column. The rightmost column shows double the size of the other two (the dotted square corresponds to the area shown in the snapshots located to the left). For the rightmost column, notice the existence of a surviving fragment at r ~ (900, −1200) au for the midplane views.

Fig. 2 (a) Mass, (b) luminosity, and (c) accretion rate of the central massive protostar, as formed in run x16. The dotted line in panel b corresponding to the stellar luminosity is calculated using the evolutionary tracks of Hosokawa & Omukai (2009). The gray line in panel c shows a reference accretion rate of 10−3 M⊙ yr−1.

4 Post-processing methods

Since fragmentation was studied without the use of sink particles, sophisticated post-processing algorithms had to be developed in order to track their properties in time, which we present in this section together with strategies for isolating the properties of the background disk.

4.1 Properties of the disk

During the fragmentation epoch, the density maps of Fig. 1 show three distinct components (background disk, spiral arms, and fragments), which we mentioned in the preceding section. In order to study the properties of the background disk, the spiral arms and fragments have to be filtered out from the data. We do this with a combination of the two following methods. First, in order to get the radial profile of a quantity q measured for the background disk in the midplane, we take the median of all values along the azimuthal direction, that is, for each discretized radial distance r_i, $$q_{\text{bg},\text{midpl}}(r_i)=\underset{k}{\text{median}}q(r_i,\theta =\pi /2,{\varphi }_k).$$(12)

Taking the median filters out most of the variability caused by the presence of fragments and spiral arms. For some radial profiles, an average of the radial profile in some interval of time is also used in order to eliminate strong but short-term variations. As fragments and spiral arms are highly dynamic, a time average filters out sudden changes and yields the long-term behavior of the background disk. With this method, one can isolate, for example, the density and temperature profiles of the primary disk, which we present and discuss in Sect. 5.1. The specific time intervals used for the average are specified in the caption of the relevant plots.

In order to determine whether or not the formed disk is Keplerian, the deviation from gravito-centrifugal equilibrium, or Keplerianity, of the background disk is calculated as the relative difference between the angular velocity Ω (taken for the background disk in the midplane) and its Keplerian value, Ω_K : $$\text{Keplerianity}=\frac{\Omega -{\Omega }_{\text{K}}}{{\Omega }_{\text{K}}},\text{where}{\Omega }_{\text{K}}(r)=\sqrt{\frac{G{M}_{\text{encl}}(r)}{r^3}}$$(13)

is the Keplerian angular velocity, and M_encl(r) is the total mass enclosed in a sphere of radius r; that is, the sum of the mass of the central massive protostar and the portion of the disk enclosed. The Keplerianity was used to define the disk radius, as is mentioned in Sect. 5.1.

The accretion rate onto the central massive protostar, as well as its mass, are calculated by integrating the incoming mass from all directions at the surface of the sink cell. One general limitation of the sink cell approach is the inability to distinguish between accretion of fragments onto the stellar surface and formation of close companions. This problem is discussed in Sect. 7, together with the implications for the values obtained with the method described here.

We remind the reader that no artificial viscosity is introduced in these simulations to compensate for unresolved physics; self-gravity instead is the process that provides angular momentum transport. This motivates us to calculate how the mass is transported by the background primary disk, that is, Ṁ (r). We take the background disk density ρ_bg(r_i, θ_j) and radial velocity v_r,bg(r_i, θ_j) and then weintegrate in a spherical shell according to $$\dot{M}(r_i)=-{\sum }_j{\rho }_{\text{bg}}(r_i,{\theta }_j)v_{\text{r},\text{bg}}(r_i,{\theta }_j)r_i^2\sin ({\theta }_j)\cdot 2\pi \cdot \Delta {\theta }_j.$$(14)

Outside of the disk, this quantity simply expresses the mass flux transported by the infalling large-scale envelope.

In order to quantify gravitational instabilities, we use the Toomre parameter Q (Toomre 1964). This parameter measures gravitational stability against small perturbations in a rotating, self-gravitating disk, and is defined as $$Q=\frac{c_{\text{s}}\Omega }{\pi G\Sigma },$$(15)

where Σ = ∫ ρ dz is the surface density and c_s is the speed of sound; Q < 1 indicates instability.

4.2 Fragment detection

In order to study the different properties of the fragments, one must first clearly define what a fragment is. One possibility is to use the fact that fragments appear as hot points in comparison with the background disk and spiral arms (see, e.g., the temperatureplot during the fragmentation epoch in Fig. 1). In comparison, the density data for the same time show spiral arms, filaments, and other structures that make the definition of a fragment more difficult. In other words, spiral arms generally have temperatures that are similar to the background disk, in sharp contrast to the fragments, allowing us to fix a criterion for fragment detection. Detecting a fragment then reduces to correctly detecting maxima in the temperaturedata for the midplane. Fragments are expected to gravitationally collapse, and therefore increase their temperature, if they have a possibility to form stars (below, we show that not all fragments can evolve further into stars).

In order to filter out the spiral arms and the disk, a Gaussian filter is first applied to the temperature data (this blurs out the spiral arms), and then the radial disk temperature profile is subtracted from the temperature data. After that, a temperature threshold is set so that the fragments are detected (~ 400 K). The parameters for such filters were manually tuned by selecting the values that produce fewer false positives and do not leavefragments undetected. This process produces “fragment candidates”. However, it should be noted that the data reveal cold high-density regions that are not considered as fragments unless they overcome the temperature-detection threshold. The fragment candidates are also checked by eye to filter out any remaining false positives that, when compared against the density data, turn out to be hot areas of spiral arms and not clump-like structures.

4.3 Fragment tracking

We are interested in calculating the properties of the fragments, and their evolution in time; therefore, we require a way to track them using the data snapshots that are output every 10 yr of evolution, which is a much smaller timescale than the duration of one orbit (≳ 500 yr at ~ 100 au). This is done as follows: for every fragment, the next predicted position on its orbit is estimated for a time t_current + Δt with Δϕ ~ v_ϕΔt. Time is then advanced by Δt, and maxima in the temperature are searched for in a region surrounding the predicted position. When a matching fragment is found in the predicted region, a connection is made and it is registered as a fragment with a unique identification number. As a warning to the reader, even though the identification number for new fragments increases with time, it cannot be used as an indication of the total number of fragments present because manual corrections to the output of our algorithm cause missing identification numbers.

4.4 Properties of the fragments

Once the trajectory of the fragments is known, we can easily track other properties, such as the central temperature or the mass of the fragments. The central temperature is the maximum value of the temperature field of a three-dimensional region that contains the fragment. The mass of the fragment M_fragm is obtained by integrating the density in a spherical region of a fixed radius (all fragments are assumed to have the same radius). A first guesstimate from the x8 density maps is R_fragm ~ 50 au.

In order to obtain an idea of the upper limit of the radius of a fragment, we can estimate the order of magnitude of this radius assuming hydrostatic equilibrium (pressure of self-gravity $P~G{M}_{\text{fragm}}^2/R_{\text{fragm}}^4$ is equal to the pressure of the ideal gas P ~⟨n⟩k_BT), using a mean number density ⟨n⟩~ 1.5 × 10¹⁷ cm⁻¹ extracted from the density data, and a temperature of 2000 K. This yields R_fragm ~ 30 au for a fragment of 1 M_⊙, and R ~ 50 au for a fragment of 2 M_⊙.

To justify our choice, we performed mass calculations with several assumed radii (the results are shown in Fig. 3). Small assumed radii give a steep gradient in the mass value, suggesting that most of the mass of the fragment has not yet been enclosed in the region of integration. At around 50 au for run x8, the curve shows a more stable value. Bigger radii would not only enclose the fragment but also parts of the disk and surrounding spiral arms. A similar analysis done for run x16 yields a radius of 40 au. These considerations mean that the mass of the fragments reported here have an uncertainty of around 20%.

5 The fragmenting accretion disk

The following sections present the results of the highest resolution simulations, and an in-depth analysis of the internal processes that govern fragmentation in the disk.

Fig. 3 Variation of the time-averaged mass of a sample of fragments as a function of the assumed radius taken for integration.

5.1 Properties of the disk

Radial profiles of several quantities related to the background disk are presented in Fig. 4, obtained using the data from run x8. The density, Keplerianity, and radial mass flow are averaged in time, as discussed in Sect. 4.

After the disk-formation epoch (from 4 kyr until the end of the simulation), the disk is approximately Keplerian, as shown by the Keplerianity profile, with a variability of ~10%. We used this fact to define the disk radius as the point at which the Keplerianity of the background disk drops below − 15%. For larger scales, the low values of Keplerianity indicate infall from the large-scale envelope, which replenishes the disk. The initial rotation profile chosen, with β_Ω = −3∕4, means that the gas is uniformly non-Keplerian, and during the disk formation epoch (0 to 4 kyr), the disk builds up until it reaches gravito-centrifugal equilibrium. In some of the panels, the disk radius is indicated as a reference with filled diamonds. The density profiles show that the background disk stays mainly between 10⁻¹⁵ and 10⁻¹³ g cm⁻³. A higher density region is observed near the edge of the disk, and corresponds to a centrifugal barrier.

As expected, the temperature profile of the collapsing cloud, including the accretion disk, increases over time, although we observe that the profile increases more slowly during the fragmentation epoch. The temperature profile is approximately proportional to r^−0.5. At the late stages of thesimulation (cf. the line for 16 kyr in Fig. 4c), there is an increase in temperature in the inner parts of the disk (≲40 au) that is probably due to stellar irradiation because the central massive protostar has just begun burning hydrogen at that time.

The radial mass flow transported by the background disk is different from the radial flow in the infalling large-scale envelope. This shows that mass from the infalling envelope is not only being transported onto the central massive protostar but is also deposited into the fragments and spiral arms, and that during the fragmentation epoch most of the mass is delivered to the massive protostar by means of the spiral arms and fragment accretion and not through the background disk.

Data from run x8 reveal a region with a seemingly unusual overdensity (more associated with the densities of spiral arms) at about 12 kyr in the region r ≲ 100 au, which is particularly evident on the surface density and radial mass flow profiles. As a result of an accretion event in run x8, the inner disk temporarily acquires a higher density “ring like” accretion structure. As a consequence, computing the background disk density as the median value in an annulus yields a higher density. Such accretion structures are also occasionally observed in the other runs, although with shorter durations. Disregarding this short-term feature, the general trends shown in Fig. 4 are also observed in the other runs.

The accretion rate onto the central massive protostar as a function of time shows great variability, because fragmentation creates small accretion events on top of having a very smooth, constant process, as described in Sects. 6.6.2 and 7.3.1. Nevertheless, the mean accretion rate in the fragmentation epoch is fairly constant, at about 10⁻³ M_⊙ yr⁻¹, although a slight increment with time in the accretion rate is observed in both panels e and f, which is due to simple acceleration of the gas infree fall.

5.2 Local versus global Toomre parameter

The Toomre maps in Fig. 1 show that the Toomre parameter is highly dependent on position and time. In Fig. 5, we present the Toomre parameter as a function of radius, averaged in time for the fragmentation epoch. The curve for Q_disk is calculated by taking the azimuthal median of all the variables and then computing Q, and corresponds to the values for the background disk. On the contrary, the curves for Q_mean and Q_min are obtained by first calculating Q(r, ϕ) and then taking the mean and the minimum, respectively. Q_min represents the regions of the disk that undergo most fragmentation, that is, spiral arms and fragments, while Q_mean may be interpreted as a value ofQ(r) obtained observationally for spatially unresolved sources (i.e., substructures such as spiral arms and fragments are not detected).

The background disk is Toomre-stable, and no fragmentation is expected to occur in these regions; however, spiral arms and fragments are subject to fragmentation. The mean value, on the other hand, shows stability. The whole disk is therefore globally stable while being locally unstable. This means that the substructures of the disk have to be resolved accurately in order to capture fragmentation; an insufficiently resolved disk may appear to be Toomre-stable while undergoing fragmentation at unresolved scales.

5.3 Spiral arm formation

In all runs, at the end of the disk formation epoch (t ~ 400 yr, see Fig. 1), a ring-shaped region in the disk becomes Toomre-unstable and develops small inhomogeneities that become two primordial spiral arms at opposite sides of the disk, that is, they arise from the m = 2 mode described in Laughlin & Rozyczka (1996) (see also, e.g., Kratter et al. 2008; Kuiper et al. 2011). The first fragments form at the outermost parts of the spiral arms.

New spiral arms are created by convergent flows, as shown in Fig. 13a; the flows left of and below the fragment are convergent and mass accumulates as spiral arms. These type of flows are frequently created by turbulent motions arising after fragment interactions.

6 Fragments

In the following sections, we examine the different physical processes that occur during the lifetime of a fragment: its formation, the formation of a secondary disk around it, the interactions with other fragments and their environment. Additionally, we present some properties of the fragments: the number of fragments present in the primary disk, their masses and their orbits, all of which were calculated with with the methods presented in Sect. 4.

Fig. 4 (a)–(e) Radial profiles of several quantities for the background disk. The intervals in profiles (a), (d), and (e) mean that the quantity has been time-averaged in that interval, excluding the endpoint. The filled diamond indicates the value at the disk radius. (f) Accretion rate into the central massive protostar as a function of time.

Fig. 5 Comparison of the Toomre parameter values computed for the background disk (Qdisk), the mean value (Qmean), and the minimum values corresponding to the effects of the fragments and spiral arms (corresponding to Qmin). These curves were calculated with run x16, and they are time-averaged for the fragmentation epoch.

6.1 Fragment formation

Fragments form out of one or more spiral arms. This two-step fragmentation process (first the disk forms spiral arms, and then they fragment) is in agreement with what Takahashi et al. (2016) reported in the context of simulations of self-gravitating protoplanetary disks. We have identified two distinct mechanisms of fragment formation in our simulations: a local collapse within a spiral arm and a formation triggered by a collision of two spiral arms. These mechanisms are present at two different scales: the primary disk, and the secondary disks that form around the fragments (we discuss the latter in more detail in Sect. 6.7).

The first process, local collapse within spiral arms, is illustrated in panels A1 and A2 of Fig. 6. When small inhomogeneities are developed within a spiral arm or when a spiral arm develops a crease, it becomes Toomre-unstable and some regions inside it begin to collapse, forming small fragments that deform the spiral arm. This process is exemplified by the creation of fragments 16, 13, and 14 of panel A3 (run x16), where the spiral arm develops inhomogeneities that grow into fragments. The other example is given in panel A4, where fragments 54 and 59 are formed by a previous crease; in the case of fragment 59, the crease was formed in a spiral arm of the secondary disk around fragment 40.

The second process, spiral arm collision, is illustrated in panels B1 and B2. When two spiral arms or regions of high density collide, theycreate perturbations that make the region Toomre-unstable and therefore, a collapse is triggered. The example in panel B3 is from run x8, and shows the formation of fragment 32, which is generated by this mechanism. Panel B4 provides an additional example from run x16, where two spiral arms, connecting fragments 32 and 68 to the inner disk, have collided and have given rise to fragment 71.

Early during the fragmentation epoch, both processes occur predominantly in the primary disk. However, towards the end of the fragmentation epoch, and as fragments gain more mass, the fragmentation inside substructures becomes dominant: the spiral arm breakdown and spiral arm collisions become more frequent in the secondary disks and yield fragmentation. This behavior was observed especially in runs x4 and x16. This hierarchical fragmentation is in principle analogous with the observations and conclusions offered in Beuther et al. (2015, 2019).

Fig. 6 Fragment creation mechanisms. The color scale for the midplane density maps is the same as the one shown in Fig. 1.

6.2 Number of fragments

The number of fragments present in run x8 as a function of time is plotted in Fig. 7. We remind the reader that the numbers next to the fragments in the plots throughout this paper are mere identification numbers. In total, 60 fragments are detected during the fragmentation epoch (4–15 kyr). From these, 22 live longer than 200 yr; these fragments will be the focus of the analysis that follows. At a given time, there are fewer than 9fragments present in the disk. A discussion about the number of surviving fragments is offered in Sect. 7.2. A peakin fragmentation is seen at around 9 kyr, although we do not consider it to be of great interest, because it is not observed in the other runs (see Fig. 20).

Lifetimes of the fragments (color code) and the radial location when they are formed (vertical axis) are shown in Fig. 7b as a function of time (horizontal axis). We see that the longest lived fragments are formed in the outer disk, and that the inner disk forms fewer fragments over time. This trend is tied to the fact that fragmentation is occurring inthe spiral arms of the primary disk. The figure was generated with the data from run x8. In run x16, as the spiral arms of the secondary disks fragment, new long-lived fragments start to form again in the inner disk.

6.3 Mass

Figures 8 and 9 (corresponding to runs x8 and x16, respectively) show the masses of the fragments obtained using the method described in Sect. 4. The masses of the fragments are of a few solar masses, and increase over time. Peaks and sudden increases are observed, and correspond to fragment mergers, which are described in Sect. 6.6. Some decreases and oscillatory behavior are caused by the action of spiral arms associated to the secondary disk, as described in Sect. 6.7. The first fragments formed tend to have lower masses than the fragments formed laterin time. This is not surprising, because the fragments form at larger radii, where larger portions of the disk become unstable, but also because the radial mass flow of the background disk increases both with the radial position as well as withtime. More massive fragments tend to have longer lifetimes as well because they are more strongly gravitationally bound and are more immune to the fragment-destruction mechanisms described in the following sections.

Two special cases occurring in Fig. 8 merit a more detailed discussion. Fragments 11 and 94 show unusually high masses;fragment 94 in particular shows high variability. Fragment 11 undergoes several mergers and gains mass just before it goes into the sink cell. The erratic behavior of fragment 94, on the other hand, is mostly numeric: it moves outwards in the simulation domain, where there is less resolution; and by effect of numerical diffusion, increases in radius, making our assumed radius of 50 au insufficient in the calculation. The mass of fragment 94, assuming a larger radius of ~ 100 au, is ~ 6.5 M_⊙ at t = 13 kyr.

6.4 Hydrostatic cores

As we discuss in Sect. 7, fragments collapse and form hydrostatic cores. We remind the reader that the masses presented in this section are calculated by assuming a fixed fragment radius of 50 au for run x8, and 40 au for run x16,which means that the enclosed region includes the hydrostatic core as well as the secondary disk. This isimportant when comparing to models of core collapse, such as the ones described in Bhandare et al. (2018): for an initial core mass of a few solar masses, they find hydrostatic first Larson cores of ~ 3 au in radius, with masses of the order of a few 10⁻² M_⊙.

In run x16, cell sizes of ≲1 au are reached at radial positions of r ≲ 100 au, barely allowing us to resolve the hydrostatic core region of fragments located there with a few grid cells, enough for an order of magnitude check with a core collapse model. We calculated the enclosed mass in the inner ≈ 3 au of fragment 12 (run x16) over its evolution. During the gravitational collapse of the fragment (Fig. 14a), the mass of the inner ≈ 3 au increased, until a value of a few 10⁻² M_⊙ was reached, consistent with the results of Bhandare et al. (2018). This high-density inner region is also shown in Figs. 13b and 14b, surrounded by the secondary disk (the figures correspond to the midplane and vertical cuts in the density field, respectively). As a reference for the reader, the yellow circle in Fig. 13b indicates the inner 3 au (roughly the size of the first core), and the yellow circle in both panels of Fig. 14 indicates the assumed size of a fragment, which was used for the calculation of its (total) mass.

Fig. 7 (a) Number of fragments with a minimum lifetime of 200 yr present in the simulation as a function of time. (b) Radial position at formation time for each fragment, color-coded with the fragment lifetime.

Fig. 8 Mass of the fragments with life span longer than 200 yr. A fixed radius of 50 au is assumed for the calculations of this plot.

6.5 Fragment dynamics: orbits

The remainder of Sect. 6 is devoted to the study of how fragments move in space, their interactions with their environment, and their substructures. First, an overview of the orbits is presented. In Fig. 10, the orbits of long-lived fragments (>1 kyr) have been plotted. During the simulated time period, only about eight fragments complete more than one orbit (both runs). The orbits of thefragments are highly influenced by their interactions with spiral arms and other fragments (fast migration, mergers, etc.; more details are given below). They also stay in the midplane, except for some orbits of run x16 that show a small inclination (≲ 2°) with respect to the midplane during certain times. The average period is ~ 1 kyr. The orbits of the fragments are highly eccentric, with an average eccentricity of ~ 0.5, which was calculated by taking the minimum and maximum radial positions of the fragments. Fragments undergo different types of interactions that cause changes in the eccentricity of their orbits, including fast inward migration due to spiral arm action and gravitational interaction with other fragments.

The orbits of long-lived fragments 40 and 32 of run x16 (Fig. 10b) develop in the middle disk,and correspond to fragments that survive at the end of the simulated time. These fragments are also more massive than short-lived ones (according to Fig. 9, their masses are ~ 1.5 M_⊙), and therefore their self-gravity provides more stability against interactions with the environment.

However, the orbits presented here do not take into account the effects of the formation of second Larson cores and the fate of fragments that enter the numerical sink cell. A more detailed discussion is presented in Sect. 7.2.

6.6 Interactions

6.6.1 Fragment–fragment interactions

Fragments interact with each other in two distinct ways: they can merge, or change orbits due to gravitational interaction and angular momentum transfer. Mergers (panel A1 of Fig. 11) occur when the fragments have a close encounter, usually in a collision orbit. As shown in the example provided in panels A2, A3, and A4 of Fig. 11, fragments 42 and 40 are in a collision orbit and merge. The masses of the fragments combine, as shown by the spike in the mass of 40 in panel A4. However, we highlight several features of a typical merger: first, the mass of fragment 42 increases slightly near the collision time. This can easily be explained by the fact that the mass integration domains start to overlap because they are spheres of a fixed radius.A second feature is the decrease in mass, with small oscillations being observed after the collision. Collisions typically do not occur head on, but with a certain impact parameter, which increases the spin angular momentum of the merged fragment. In the frame of reference of the fragment, the increased centrifugal force favors the development of secondary spiral arms, which can transport mass outwards, as explained in Sect. 6.7. During the fragmentation epoch in run x8, we recorded about 13 mergers in total. This kind of interaction can only be captured by simulating and resolving the full hydrodynamics of a fragment, and not by a sub-grid particle model.

Changes in the orbit due to gravitational interaction of fragments that are unconnected by a spiral arm are somewhat more complex, less frequent, and more difficult to determine. However, a general picture is provided by panel B1of Fig. 11: two fragments approach and their gravitational interaction slows down fragment 1, therefore moving it to a lower orbit, and accelerates fragment 2, moving it to a higher orbit. Panel B4 provides an example of this kind of interaction: just after the orbits of fragments 16 and 29 cross each other, fragment 16slows down and enters into a collision orbit with the central massive protostar, while fragment 29 gets accelerated into a higher orbit. However, the moments leading to the orbit crossing are another example of the effects of a gravitational interaction: in panels B2 and B3, we see how fragment 29, originally in a low orbit, migrates towards fragment 16 and vice versa, until they cross each other and the scenario illustrated as t₀ in panel B1 is reached. The example provided here comes from the data of run x4, but we also observed a similar behavior in run x8, but without the orbit crossing. Some of these interactions have also been observed previously in 2D simulations in the context of planet formation (see, e.g., Zhu et al. 2012). The outcome of these interactions (a change in the orbit or a merger) depends on the approach velocities and the mass ratio of the interacting fragments. In the example given in panels B2–4, the masses were comparable.

Fig. 9 Masses of the fragments with life span longer than 200 yr, for run x16. A fixed radius of 40 au is assumed for the calculations of this plot.

Fig. 10 (a) Orbits of fragments with a life longer than 1 kyr during the time interval [6, 10] kyr, for run x8. (b) Orbits of all the fragments with a life longer than 1 kyr during all of the simulated time for run x16. The fragment identification number is located at the starting point of the orbit, and the color code is the same as that used in Figs. 8 and 9.

6.6.2 Fragment interactions with spiral arms and the massive protostar

The gravitational torques exerted by the spiral arms can cause migration, as shown by panels A1–4 of Fig. 12. In this case,the eccentric orbit of fragment 0 reduces its periastron due to the action of the connecting spiral arm.

Accretion onto the central massive protostar is not a smooth process, as depicted in the accretion rate shown in Fig. 4f. There are many discrete accretion events caused by infall of matter through a spiral arm or the complete accretion of a fragment. The transformation of gravitational energy into radiation is made in the form of accretion-driven bursts, that is, sudden increases in the luminosity of the forming massive star (see Fig. 2b). This phenomenon was reported in theoretical studies in the context of massive star formation by Meyer et al. (2017, 2019), where a system of magnitudes was also developed to describe them.

Panels B1–4 of Fig. 12 show two examples of accretion bursts: both the infall of matter through a spiral arm, as shown in panels B1–2, and the complete accretion of a fragment (panels B3–4) produce accretion bursts. However, the accretion of a fragment produces a much sharper and brighter peak in luminosity compared to accretion through a spiral arm. Accretion bursts are also accompanied by overall increases in the temperature of the disk, which we term temperature flashes, which typically last ~ 30 yr and heat up the fragments and spiral arms, as we discussed in the previous section.

Nevertheless, we must discuss the effects of the inner boundary conditions (size of the central sink cell) on these accretion bursts. As matter is only allowed to cross the sink cell inwards, we are not taking into account how outflows affect the accreted mass. Some fragments that the simulation shows as accreted into the sink cell could also only be in an elliptical orbit, as discussed in Sect. 7, especially if they have already undergone second collapse, meaning that they are less susceptible to shearing and subsequent accretion. These considerations could make strong luminosity accretion bursts less frequent and intense than presented here.

Fig. 11 Fragment–fragment interactions. Panels A3 and B4: the fragment ID is shown at the starting point of the orbit. Panel A3: the time window shown in the orbits is 1.3 kyr; panel B4: this is 1.2 kyr.

Fig. 12 Interactions of the fragments with the spiral arms, and with the forming protostar modeled as a central sink cell. The orbit in panel A4 is plotted for a window of 2 kyr, with the fragment ID marking the starting point.

Fig. 13 Secondary disk with its spiral arms around (a) fragment 11 of run x8, and (b) fragment 12 of run x16. The background map represents density in the midplane (using the same color scale as the rest of the figures), and the arrows show the comoving velocity field. The central massive protostar is located beyond the bottom of the plotted areas in both cases. The yellow circle in (b) has a radius of 3 au.

6.7 Secondary disks

Fragments frequently develop a secondary disk with spiral arms, which is embedded in the primary disk, as shown in Fig. 13. The figures show the comoving velocity field, v − v_ϕ,bge_ϕ, where v_ϕ,bg is the angular velocity profile of the primary disk. In Fig. 13a, we observe a typical secondary disk, with the densest and hottest region (fragment) in its center. We also observe how converging flows form spiral arms. Converging flows form frequently when shearing and turbulent motion (due to activity of fragments and other spiral arms)produce a region with a net outward flow that encounters the inward radial flow from global mass transport.

Spiral arms in secondary disks can also be drifted off the fragment, a phenomenon exemplified in Fig. 13b: the spiral arm encounters a decelerating inward flow that makes it drift apart. As a result, the calculated mass plot (Fig. 9) shows a decrease because the mass of the spiral arm is lost from the fixed integration domain used. If the spiral arms are bounded to the fragment, the plotted mass of the fragment only shows an oscillation, but frequently the spiral arm gets absorbed by other spiral arms, provoking a real mass loss in the fragment with a consequent temperature variation over time. This mechanism provides a way to transfer some matter from fragments back to the disk or the central massive protostar. This process can also trigger fragmentation, as was the case of the spiral arm that gives rise to fragment 59 shown in Fig. 6A4.

Figure 14 shows the vertical structure of a fragment using the highest resolution dataset. When the fragment is forming (Fig. 14a), it has an approximately spherical shape. However, Fig. 14b shows the vertical structure of the secondary disk shown in Fig. 13b: the fragment indeed gets flattened down, with the highest density in the center.

6.8 Comparison to observations

6.8.1 G11.92–0.61 MM 1

Observations of the massive young stellar object G11.92–0.61 MM 1 reported in Ilee et al. (2018) show a fragmented Keplerian disk around a proto-O star. The main object, MM1a, is reported to be ~ 40 M_⊙ (with between 2.2 and 5.8 M_⊙ attributed to the disk, and the rest to the protostar), and the fragment MM1b is reported to be ~ 0.6 M_⊙, located at around ~2000 au from MM1a.

This seems to be especially compatible with the general results in run x8, but later in time (the disk reported in the observationsis bigger, and the central star is more massive). However, there are several warnings that should be taken into account. From the results in Ilee et al. (2018), the mass ratio of the disk and the primary is between 0.055 and 0.145. In our simulations, the mass of the central massive protostar at ~16 kyr is ~ 20 M_⊙ (cf. Fig. 19). We calculated the mass of the disk including the fragments and spiral arms (integration of density in a cylinder of 1500 au in radius and 20 au in height), and the mass of the background disk. The mass of the disk including substructures is comparable to the mass of the primary, as expected from a fragmenting disk; however, if substructures are excluded, the mass of the disk increases more slowly. After ~ 16 kyr, the mass of the background disk is ~ 5 M_⊙, which means a mass ratio of ~ 0.25. In addition, the surviving fragment (94) in run x8 is located at ~1500 au from the primary, although at the end of the simulation it is moving outwards. The mass of fragment 94 is of the order of a few solar masses, although this value includes the mass of the secondary disk.

Given that the masses of most of the fragments produced in our simulations are of the order of 1 M_⊙, and given the position and size of the disk both reported here and in the observations, we think that the scenario described in Ilee et al. (2018) could plausibly be obtained with a setup similar to ours.

Fig. 14 Vertical structure of fragment 12 and its close surroundings (run x16). The yellow circle has a radius of 40 au. The density uses the same color scale as the rest of the figures.

6.8.2 Accretion burst event in G358.93-0.03

Disk substructures associated to an accretion burst were observed and reported in Chen et al. (2020) for the high-mass young stellar object G358.93-0.03. The flaring event was reported in Sugiyama et al. (2019). Chen et al. (2020) performed a kinematic model that describes the accretion flow as occurring along two spiral arms, although they also acknowledge the compatibility of their observations with the accretion of a fragment. The results of our simulations support the second scenario, as do the results in Meyer et al. (2017, 2018). The processes described in Sects. 6.5, 6.6.1, and 6.6.2, namely the gravitational interactions between fragments (as exemplified in panel B4 of Fig. 11), and the gravitational torques exerted by the spiral arms can cause the fragment to lose angular momentum, leading to accretion by the primary, and thus an accretion burst.

Fig. 15 Temperature of the fragments with lifespan longer than 200 yr, for run x8. The gray box indicates the hydrogen dissociation limit.

7 Companion formation

As fragments contract by their self-gravity, their temperatures increase, which in turn causes an increase in their internal pressure, halting the collapse (classical Kelvin-Helmholtz contraction). However, fragments whose central temperature goes beyond ≈ 2000 K start to experience hydrogen dissociation, which means that gravitational energy is not converted into thermal kinetic energy, but rather is used to dissociate hydrogen molecules, therefore allowing for further (second) collapse, until a second Larson core is formed.

The exact temperature for hydrogen dissociation is density-dependent; we verified that for the central densities of the fragments, the temperature is between 1700 K (10% dissociation) and 2300 K (90% dissociation) (D’Angelo & Bodenheimer 2013). According to Bhandare et al. (2018), second Larson cores formed from reservoirs of a few solar masses have radii of the order of a few solar radii. Based on the free-fall timescale of a fragment of a few solar masses, we use the conservative estimate of 80 yr for the durationof the second collapse. After the second core is formed, contraction continues to drive the temperature up, continuing the evolution of the protostar.

The simulations presented here do not include hydrogen dissociation. However, by tracking the central temperature of a fragment, one can determine if it should undergo second collapse and form a second core. This allows us to hypothesize the ultimate fate of the fragments if more realistic physics were considered.

7.1 Central temperature of the fragments

Figures 15 and 16 show the tracked central temperature of the fragments as a function of time, for runs x8 and x16, respectively. The shadowed box shows the threshold for hydrogen dissociation, which, as mentioned above, is not included in our simulations. This is why we see high temperatures and fragments crossing the threshold both ways over time. In both figures, we see some high and short-lived spikes in the central temperature, which are the result of rapid mass increase, for example due to mergers or accretion onto the fragment. The small periodic variations of the temperature are caused by the interactions between the fragments and their environment, including the development of secondary spiral arms. The mechanism of secondary spiral arm drift explained in Sect. 6.7, which produces a lowering in the mass of the fragment, also produces a lowering in its temperature (due to decompression). The variation of the observed central temperature of the fragments with resolution is discussed in Sect. 8.2.

Longer lived fragments tend, in general, to have higher temperatures of at least ~ 1000 K. Fragments formed later in time tend to also reach higher central temperatures, because they have more mass (cf. Sect. 6.3).

Fragments that cross the hydrogen dissociation temperature threshold undergo second collapse and become second cores. As a reference, we count the number of fragments that have a temperature higher than ~ 2000 K for a time longer than ~80 yr, which is ourestimate of the time needed to form a second core. In run x8, four second core objects are formed, and in run x16, ten second core objects are formed.

7.2 Fate of the fragments

Only after the second collapse has taken place would a fragment reduce its radius from the order of a few astronomical units to a few stellar radii. However, merely from the data provided by the simulations presented here, it is impossible to know whether or not the temperature continues to rise after the second collapse and whether or not the second cores become actual companion stars. Nevertheless, one point to take into account is that many fragments develop secondary disks, which means that after the second collapse, more mass can be delivered to the second core and potentially allow a temperature increase and continuation of the evolutionary process.

Fragments of sizes of a few solar radii would be unresolvable with our current simulation setup. However, since their compactness would also provide stability, it would be safe to add a sub-grid particle model in order to study their behavior. However, before the second collapse happens, fragments are mere hydrostatic cores with strong interactions with the disk, spiral arms, and other fragments, and some of these interactions are strong enough to destroy them. The following sections describe these mechanisms, and consider the effects that our simulation setup with a central sink cell have on the number of companions that are formed according to the simulations.

7.2.1 Fragment destruction mechanisms

Fragments that move in an eccentric orbit near the central primary protostar with high speeds experience shearing, which typically causes their destruction, as illustrated in panel A1 of Fig. 17, and exemplified by panel A2. A similar process can be found in simulations of fragmentation in protoplanetary disks, as described in Lichtenberg & Schleicher (2015). It is also common in this scenario that some matter of the fragment gets accreted by the central protostar, sometimes forming a “ring-like” accretion structure like the one discussed in Sect. 5. However, in the particular example of panel A2, matter is not accreted. The sheared material expands and becomes absorbed by spiral arms, providing converging flows needed for further fragmentation according to the mechanisms discussed above.

Spiral arms can also cause mass loss from the fragment by transporting matter along them (Fig. 17B1), eventually destroying the fragment. As an example, fragment 43 of run x8, shown in panel B2, is drained by a spiral arm that is being stretched, and delivers the material to fragments 32 and 40.

Another fragment-destruction mechanism occurs as a consequence of temperature flashes associated with accretion bursts (described in Sect. 6.6.2). Fragments and spiral arms are heated up by the temperature flash. This triggers thermal expansion, which in turn lowers the density and central temperature of the fragment, compared with the values before the flash. If the fragment already has a low temperature and density, this thermal expansion and cooling can take the fragment below the detection threshold. After manually following the remaining matter, we see that it gets partially dissolved into the disk, absorbed into a spiral arm, absorbed by other fragments, or is sheared apart. An example of this behavior is provided by panels C1 and C2 of Fig. 17.

These destruction mechanisms severely affect the lifetime of fragments if they are in their hydrostatic core phase, and show that sub-grid sink particle models may overestimate the number of formed companions. Only by resolving the substructure of the fragments do these interactions become apparent.

Fig. 16 Temperature of the fragments with lifespan longer than 200 yr, for run x16. The gray box indicates the hydrogen dissociation limit.

Fig. 17 Fragment-destruction mechanisms.

7.2.2 Fragment dynamics after the second collapse

Once a second core is formed, the question arises as to how such an object will interact with the environment, and therefore, how the fragment destruction and merger mechanisms discussed above will impact the fate of a fragment. Although second cores would experience the same gravitational interactions as a hydrostatic core, they would not experience the same pressure gradient from the rest of the gas due to their compactness. A second core would have less probability of colliding and merging with another second core because of the reduction in the collisional cross section. A merger between a first core and a second core is a possible scenario: the first core would likely be sheared apart, and would form a denser secondary accretion disk around the second core, ultimately being accreted by the latter.

7.2.3 Formation of spectroscopic multiples

As mentioned in Sect. 4, one disadvantage of using a central sink cell is the inability to distinguish between accretion and the formation of close companions. This is especially true for fragments that have undergone second collapse, and therefore donot suffer the destruction mechanisms mentioned above. Once a fragment enters the numerical sink cell (r = 30 au), its mass is counted as being accreted by the central massive protostar, when in reality, if the fragment is a second core, it might be in an orbit instead of a merger with the primary.

Migration and gravitational interactions with other fragments are some of the responsible mechanisms for getting a fragment into an orbit closer to 30 au (the radius of the sink cell). In order to study the possible fate of a second core, we integrated the orbit of fragment 2 of run x16 (Fig. 18), which undergoes second collapse shortly before it enters the sink cell. In the integration, we considered the gravitational force of the central massive protostar and the gravitational force of the background disk, but no other interactions.Fragment 2 was originally dragged into the sink cell by gravitational torques arising from spiral arm interaction. The predicted integrated orbit is elliptical, and it gets smaller in time because of the increase in mass of the central massive protostar. Some “jumps” can also be seen in the orbit, which are caused every time the central massive protostar accretes (non-second core) fragments, that is, when its mass suddenly increases by a discrete amount. After 3 kyr of integration, the orbit of the fragment has a periastron of ~15 au. At the time, the simulation is at t = 7.9 kyr, and the mass of the central protostar is ~4 M_⊙. As time progresses and the central protostar becomes more massive, the orbit of the fragment should become smaller and smaller, thus providing a mechanism for the formation of spectroscopic multiples.

Fig. 18 Integrated orbit of fragment 2 of run x16 after it enters the numerical sink cell. The orange dot indicates the initial position of the fragment, and the total time integrated was 3 kyr. The sink cell (dashed circle) has a radius of 30 au.

7.3 Implications

7.3.1 Accretion bursts

In run x8, fragment 11 reaches the hydrogen dissociation temperature at 10.22 kyr and would undergo second collapse if the simulation had implemented it. At the time, the mass of the fragment is 1.8 M_⊙. At 11.47 kyr, fragment 11 then merges with fragment 71, which greatly increases the temperature (cf. Fig. 15), and its mass briefly reaches 5 M_⊙. This would be an example of a second core merging with a hydrostatic core, and we would expect, in a more realistic setup, that fragment 11 would at least partially accrete the material from fragment 71 through the secondary accretion disk. However, the interaction between fragments 11 and 71 alters the orbit of fragment 11 and it ultimately reaches the sink cell. On the other hand, a second core entering the sink cell might produce a companion instead of an accretion event. For the plots of run x8 presented in this paper, this means that, at t ~ 11.5 kyr, (i) the high peak in the accretion rate (Fig. 4f), the ≳10⁶ L_⊙ accretion burst shown in Fig. 19b, and the bump in the mass of the central massive protostar (Fig. 19a) might not take place; and (ii) the subsequent inner disk overdensity that arises in the radial profiles for run x8 (Figs. 4a, b and e) during the interval [12, 16[ kyr as a result of theshearing and accretion of fragment 11 might not be there.

According to the preceding discussion, very high luminosity bursts are less likely to occur than what is shown in Fig. 2b, because they require the accretion of a significantly large mass in a short time. Fragments that have the required mass have most likely undergone a second collapse, and, according to the integrated orbits we presented above, they will likely form close companions instead of being accreted. The possibility that a second core be accreted is not excluded, but it should be rare based on collisional cross section considerations.

Fig. 19 (a) Mass and (b) luminosity of the central massive protostar, as formed in run x8.

7.3.2 Number of companions

Taking into account the preceding elements, we summarize here the results of our simulations regarding companion formation. There are four second core objects produced in run x8. From them, one (fragment 94) survives in the outer disk (r ~ 1400 au at the end of the simulation). Fragments 0 and 11 form close companions, and fragment 37 merges with another fragment (so its fate is unclear from our simulations).

In run x16, ten second core objects are produced. Fragments 0, 12, and 13 form close companions, and fragments 40, 32, and 68 survive the simulation in the middle region of the disk (see discussion in Sect. 6.5). The rest merge with other fragments, and therefore, again, their fate is uncertain from our simulation.

By analyzing which of the fragments that get close to the central massive protostar produce accretion bursts indeed, and which ones satisfy the conditions for further stellar evolution, we find that the number of possible close companions produced in our simulations is low. This result is consistent with observations of more evolved systems, such as Sana et al. (2012), Kobulnicky et al. (2014), and Dunstall et al. (2015), where spectroscopic binaries or multiples are present in large fractions, but they also do not observe tens or hundreds of spectroscopic companions to a central massive star. Nevertheless, our results show that companions produced by disk fragmentation can also exist at distances of the order of 1000 au astronomical units (although inward and outward migration is still possible in the long term). We also note in this comparison that the masses of the fragments given in the preceding sections are not necessarily an indication of the masses of the future companions, as accretion is still ongoing at the end of our simulations.

Fig. 20 Convergence of different quantities with resolution. The midplane density data were taken at t = 7 500 yr for all the maps in the row. The Toomre Q was time-averaged during the fragmentation epoch.

8 Resolution effects

After describing the main results, the focus of the following section is to describe how these change with resolution by comparing five runs with different spatial resolution.

8.1 Convergence of fragmentation

Here, we refer to the panels in Fig. 20. The midplane density maps (Fig. 20) show that only spiral arms are produced in runs x1 and x2, but no fragments. By no fragments we mean that our conditions for fragment detection are not met, and that the shape of the substructures formed (regions of higher density than the background disk) is difficult to be recognized as a fragment because of resolution. Runs x4, x8, and x16 produce fragments and spiral arms. Runs x1, x2, and x4 show thicker spiral arms than runs x8 and x16, probably as a result of numerical diffusion. Substructure formation in runs x1 and x2 continues until the end of the simulated time, that is, the evolution of the system does not drastically change over time. In contrast, a quiescent epoch with no further fragmentation is observed in runs x4 and x8. Due to computational costs, data for run x16 are not available beyond t = 13.5 kyr, which prevents us from confirming the existence of a quiescent epoch in that run as well. In summary, only the runs with sufficiently high spatial resolution exhibit an end to the fragmentation epoch.

The Toomre Q parameter shows a similar behavior for all resolutions. We see in Fig. 5 that the value of Q for fragments and spiral arms, corresponding to the Q_min, reveals the susceptibility of these areas to fragmentation, but the background disk and mean values (expected to be dominant in underresolved observations) show the system as Toomre-stable.

The number of fragments present at a given time in the simulation that live longer than 200 yr is also similar, except for the peak of eight fragments seen in simulation x8, which is discussed above. The average number of fragments at a given time during the fragmentation epoch is around 2.2. The tendency of longer lived fragments to be created at larger radial positions with time is also observed in the other simulations, although not shown here.

Fig. 21 Fragment statistics for the fragment-producing runs. The charts show the fate of the fragments with a life longer than 200 yr according to the simulation, that is without taking into account the effects of the inner boundary and hydrogen dissociation, as discussed in Sect. 7.2. The estimated number of these fragments that undergo second collapse is given in the right-hand column.

8.2 Convergence in the properties of the fragments

The fate of the fragments without taking into account the effects of the inner boundary and hydrogen dissociation we discussed in Sect. 7.2 are shown in Fig. 21. Mergers and shear/accretion into the sink cell dominate the mechanisms of fragment destruction. Spiral arm matter transport and dissolutions are far more rare.

Comparing Figs. 15 and 16, the central temperature of fragments gets higher with resolution. From this, the number of fragments that get to the second collapse phase increases from run x8 to run 16, but remains in the same order of magnitude. In Bhandare et al. (2018) (Fig. 2c), the temperature profile of a first core shows that the maximum temperature is reached at a radius of ~ 1 au. Run x16 has a resolution of ≲1 au at r ≲ 100 au, which might indicate that the values of the central temperature for fragments in the inner disk are close to convergence.

The statistics of fragment interactions show that independently of resolution, fragments in a hydrostatic core stage are fragile, and it is not safe to include them in a particle sub-grid model if they do not reach second core status, for which the correct resolution is needed. As already discussed in Sect. 6.3, the masses of the fragments for runs x8 and x16 are similar, and although not shown here, the masses of the fragments produced in run x4 also have consistent values.

Fragments in the outer disk, as exemplified by fragment 94 of run x8, tend to have higher masses, and therefore higher probabilities of reaching the temperature required for a second collapse. However, the spatial resolution of the numerical grid is lower in the outer disk (for run x8, the cell size is ~ 20 au at ~ 1000 au); and although we show in Sect. 9.2 that, at that radial distance, this is enough to resolve the Jeans length and therefore the formation of fragments, it is not enough to resolve the substructure of the fragments. Therefore, fragments located at r ≳ 1000 au have radii that are artificially higher due to numerical diffusion, even when not considering the physics of the second collapse. A larger radius makes the fragment less gravitationally bound, and filaments and secondary spiral arms are developed between the fragment and the central massive protostar. These structures fragment and the fragments migrate rapidly inwards and produce accretion bursts in the late stages of runs x4 and x8 during the quiescent epoch; but we believe this effect to be mainly numerical. Thanks to better spatial resolution, we can see that fragments produced in run x16 have more consistent sizes across the disk.

Fig. 22 Convergence of the properties of the formed protostar.

8.3 Convergence of the formed massive protostar

During the simulation, several properties of the central massive protostar are calculated under several assumptions. Figure 22 shows the mass and total luminosity. We remind the reader that the mass of the central protostar is simply the mass of the sink cell, and the total stellar luminosity is the sum of the accretion luminosity (total conversion of gravitational potential energy of the accreted material into radiation) and the stellar luminosity calculated following the evolutionary tracks of Hosokawa & Omukai (2009) (which depend on the mass of the protostar and the accretion rate).

Figure 22a shows only qualitative convergence in the protostellar mass after ~ 11 kyr, because it depends heavily on how the material is accreted in each simulation, and this in turn depends on the fragmentation that is developedin time. Figure 22b shows that the total luminosity curves are similar for all simulations: all feature some accretion bursts and a general increase in time due to the contribution of stellar evolution. The plots for runs x4 and x8 also show a general reduction in accretion bursts at the end of the simulated time (quiescent epoch), except for some bursts caused by the effect discussed in Sect. 8.2, that is, triggered by fragments in the outer disk. There is no quiescence observed for runs x1 and x2. We also find that the central massive protostar starts burning hydrogen at ~ 12 kyr for runs x16 and x8, ~18 kyr for run x4, and ~16 kyr for runs x1 and x2.

9 Comparison with previous studies

9.1 Overview

Several studies on massive star formation, including cloud and disk fragmentation, have been performed over the years. We summarize some of them in this section and compare their outcomes and main features in the following sections.

First, we discuss some studies on cloud fragmentation. Girichidis et al. (2012) performed several gravito-hydrodynamical simulations to study a collapsing cloud under the fragmentation-induced starvation scenario without radiation transport, an AMR grid of a minimum size of 13 au, and sink particles. In these simulations, various initial density profiles were used, together with a supersonically turbulent initial velocity profile but no rotation. These latter authors find between 161 and 429 sink particles in total; some of their simulations show filamentary accretion, and some the formation of an accretion disk.

Peters et al. (2010, 2011) used the code FLASH to study the effects of magnetic fields and ionizing radiation in a large-scale collapsing cloud of 1000 M_⊙, neglecting the thermodynamics of dust re-emission, that is, continuum radiation transport was not taken into account. Both studies used AMR with a minimum cell size of 98 au, and sink particles. However, according to our results, these sink particles could not be comparable to our fragments, since their accretion radius of ~ 400 au is at least one order of magnitude bigger than the size of a fragment with the potential of forming a companion, and more comparable to the size of the whole primary accretion disk, that is ~ 400 au in radius in the middle of the fragmentation epoch. With no magnetic field, around 25 sink particles were created; the presence of magnetic field reduced fragmentation by around a factor of two, and generated a more massive protostar due to magnetic braking. Similar results were obtained by Hennebelle et al. (2011), where the RAMSES code was used to study a cloud of initial mass 100 M_⊙ and radius of 1.35 pc. These latter authors used AMR with a minimum cell size of 8 and 2 au for their low- and high-resolution runs, respectively, and no sink particles, but used the barotropic equation of state. The number of fragments obtained is of the order of 50 without magnetic field.

The following studies also considered disk formation and fragmentation. In Krumholz et al. (2007, 2009), three-dimensional gravito-hydrodynamical simulations were performed, with a gray flux-limited diffusion approximation for stellar radiation as well as radiation by the dust. The first study dealt with the collapse of a core in initial supersonic shock-dominated turbulence, while the latter consiered solid-body rotation without any initial turbulence. Both studies used sink particles following the Jeans criterion only. In the turbulent case for an initial mass of 200 M_⊙ (Krumholz et al. 2007), the spherical cloud forms filaments that feed two clumps that ultimately become one massive protostar with an accretion disk that fragments into spiral arms, but no further fragmentation; however, three sink particles are independently formed at around 3000 au from the primary. With regards to fragmentation, in Krumholz et al. (2009), the disk forms earlier, and ultimately the system develops into a binary with a separation of ~ 1500 au. Both simulations used the code Orion, with an AMR grid with a maximum resolution of 10 au.

Klassen et al. (2016) used the FLASH code with an improved treatment of radiation transport, and an AMR grid with a minimum cell size of 10 au. The cloud was initially in solid-body rotation, and masses of 30, 100, and 200 M_⊙ were considered. They used a stricter sink particle algorithm that, in addition to the Jeans criterion, checks for a convergent flow, a gravitational potential minimum, and a negative total energy. The use of this criterion with the resolution considered leads to no fragmentation at all; only the formation of spiral arms was observed.

Rosen et al. (2016) repeated the simulations for both the initially laminar and highly supersonic turbulent cases, including a hybrid radiation transfer method that properly treats the multi-frequency stellar irradiation and gray flux-limited diffusion thermal re-emission (Rosen et al. 2017), but the maximum resolution of their AMR grid was 20 au. Their sink particle algorithm detects 29 companions, 16 with masses >0.1 M_⊙ for the laminar case, and 3 companions for the turbulent case. In both cases, an accretion disk with spiral arms is formed around the massive protostar, but in the turbulent case the disk becomes eccentric. Contrary to Krumholz et al. (2009), the system does not form a binary (in part due to different sink particle merging criteria). Initially, in the turbulent case the massive protostar is fed primarily by filaments, but the accretion rate is not significantly different from that in the laminar case. A similar setup was used by Rosen et al. (2019) to study the role of turbulence in fragmenting cores with initially virial and subvirial initial conditions. These latter authors find that virialized cores undergo significantly more turbulent fragmentation at early times than subvirial cores. In both cases, a fragmenting accretion disk was formed.

Meyer et al. (2017, 2018) used Pluto with a similar setup to ours, that is, a fixed spherical grid with the radial coordinate increasing logarithmically. They also used the same radiation transport scheme, but lower spatial resolution. In Meyer et al. (2018), several initial angular velocity profiles were examined for fragmentation, specifically β_Ω = 0, − 0.35, and − 0.75. These latter authors found formation of spiral arms followed by fragmentation; the highest number of fragments was obtained for β_Ω = −0.75. It was proposed that fragmentation might explain the high spectroscopic binary fraction of massive stars. Accretion bursts were observed in Meyer et al. (2017, 2018, 2019). Meyer et al. (2019) also studied the effects of disk wobbling during the fragmentation epoch.

9.2 Resolution of the grid and the use of sink particles

In order to check whether or not a simulation is resolving the correct scales associated with the physical phenomena described, we computed the pressure scale height of the disk and the Jeans length using the data from run x16, averaged over the time period [6, 8] kyr, which corresponds to the fragmentation epoch. The pressure scale height gives an idea of the scale of the vertical structure of the disk, and it is defined in the midplane as $$\frac{H_{\text{disk}}}{r}=\frac{c_{\text{s}}}{v_{\varphi }},$$(16)

where $c_{\text{s}}^2=\gamma \partial P/\partial \rho$ is the sound speed. The Jeans length λ_J indicates the size of a region that becomes gravitationally unstable at a certain density ρ and temperature, and is calculated as $${\lambda }_{\text{J}}=c_{\text{s}}\sqrt{\frac{\pi }{G\rho }}.$$(17)

In order to adequately study fragmentation numerically, the Jeans length has to be resolved properly. In order to study a Jeans length radial profile, we separate it into two values: the Jeans length of the background disk ${\lambda }_{\text{J}}^{\text{disk}}$, calculated with the median values (as described in Sect. 4), and the minimum Jeans length, ${\lambda }_{\text{J}}^{\text{min}}$, which corresponds to the value of the most fragmentation-prone areas of the disk, that is, the fragments and spiral arms.

Figure 23 contains the radial profiles of these quantities, which were calculated using the data from run x16 and time-averaged during the period between 6 and 8 kyr. The disk radius is also shown for reference. Additionally, the cell size of our simulations, Δx(r), is shown for each run (gray lines).

The pressure scale height scales approximately linearly with distance, giving an ansatz on the required cell size needed to resolve the vertical structure of the disk. The cell size of our logarithmic grid also scales linearly with distance. As a result of the simulation, we obtain H(r) ≈ 0.12r^1.07 on the disk region, and a comparison with the different runs allows us to conclude that H is resolved by all of them, but only partially by x1.

The Jeans length decreases with increasing density, which means that it increases with distance. Our choice of the coordinate grid then makes possible a better resolution of the Jeans length while saving computational power. The Jeans length for the background disk is resolved by all runs, and indeed we observe some sort of structure formation (spiral arms and/or fragments) in all runs. However, when observing the curve for the minimum value of the Jeans length, it is clear that only runs x8 and x16 are able to resolve fragmentation adequately, and run x4 only partially resolves fragmentation, which coincides with the results found in the convergence study presented in Sect. 8.

In order to compare the results from other studies, we establish approximate equivalences between the different grids used. In the case of the fixed grids used in Meyer et al. (2017, 2018, 2019): (a) Meyer et al. (2018) use approximately the same grid structure as our run x4 (except in their “Run 1-HR”, where the grid used is almost our run x8); (b) Meyer et al. (2017, 2019) use a grid comparable with our run x4.

A direct comparison of our fixed grid with AMR implementations is not trivial, because they are dependent on the selected refinement criteria and the particular evolution of the system. In order to present a panorama of the resolution differences, we take the criteria given in for example Rosen et al. (2016), that is, a Jeans length should be resolved by at least eight cells. The resulting profiles are shown in green in Fig. 23, where if eight Jeans lengths cannot be resolved, the maximum refinement is plotted.

We now examine the disk region. According to our results, the minimum Jeans length is only marginally resolved by an AMR simulation that has a minimum cell size of 5 au. Even in that case, our simulations provide much more detail, especially in the region 100 au ≲ r ≲ 600, where according to Fig. 7b we observe the biggest number of long-lived fragments being created. This choice has allowed us, as we have shown in the results, to study even the interactions between fragments and their internal structure, and not merely to resolve where a fragment should be created by the Jeans criterion.

Based on these considerations, we can establish that an AMR grid with a minimum cell size of 5 au provides a resolution between those of runs x16 and x8 for the outer disk, but with a resolution similar to or lower than that of run x4 for the inner disk (r ≲ 100 au). A minimum cell size of 10 au corresponds to a resolution that is between those of runs x8 and x4 for the outer disk, but is equivalent to that of run x2 for the inner disk. A minimum cell size of 20 au provides a resolution that is between those of runs x4 and x2 for the outer disk, but that is similar to or lower than that of run x1 for the inner disk.

Studies that use a sink particle algorithm typically use the Jeans density as the main criterion for particle creation. This is the case in, for example, Krumholz et al. (2009, 2007), Rosen et al. (2016), Peters et al. (2010, 2011). As mentioned above, a self-consistent treatment of the fragmentation process – in contrast to a sink particle algorithm – provides more certainty in the number of fragments formed, and, as we show here, in the number of fragments that survive the fragmentation process and can become companions. We also show that fragments, while in the hydrostatic core phase, undergo several interactions that can lead to their destruction, and they need to reach a second core status in order to be safely replaced with a sink particle. The use of sink particles for hydrostatic cores additionally misses the formation of the smaller secondary disks, and the fact that their associated spiral arms can occasionally give rise to new fragments.

The fact that AMR grids with minimum cell sizes of 20 au and 10 au do not adequately resolve the Jeans length, in addition to the fact that we do not observe fragmentation in runs x1 and x2, suggests that sink particles created in simulations that use such grids, although representing regions of high density, do not (necessarily) represent physical local collapse. This means that companions produced in this way are due to purely numerical effects, and therefore we do not recommend the use of a sink particle algorithm in multiplicity studies if not used in conjunction with a grid of adequate resolution.

Fig. 23 Comparison of midplane grid cell sizes in our fixed grid and other studies that use AMR. The different relevant scales of thesystem are shown as dashed lines; these were calculated as averages for the time period [6,8] kyr, and the corresponding disk radius is shown in the gray box. The number following the label for the AMR curves indicates the minimum cell size in astronomical units.

9.3 Additional physical effects

Contrary to our setup, the simulations from Krumholz et al. (2009), Rosen et al. (2016) (laminar), and Klassen et al. (2016) do not restrict the motion of the central massive protostar via a fixed sink cell. However, these studies show qualitatively that if the system forms a disk, the central massive protostar does not excessively move from the center. This was studied in more detail and confirmed in Meyer et al. (2019), where the authors implemented disk wobbling and found no dramatic differences when treating the massive protostar as fixed during the fragmentation period.

With regards to the initial density profile, studies that consider a constant density profile (e.g., Peters et al. 2010, 2011; Girichidis et al. 2012) find that the cloud forms many disperse fragments, while clouds with power-law density profiles (e.g., Girichidis et al. 2012; Rosen et al. 2016; Meyer et al. 2018, and this work) are dominated by a central object, with increasingly less fragmentation as the profile steepens. Observational studies (e.g., Beuther et al. 2002; van der Tak & Menten 2005; Palau et al. 2014) typically find density slopes β_ρ of between − 1.5 and − 2.6.

From the studies considered in Sect. 9.1, Girichidis et al. (2012), Rosen et al. (2016), Hennebelle et al. (2011), and Krumholz et al. (2007) performed simulations with supersonic initial turbulence, while Krumholz et al. (2009), Rosen et al. (2016), Peters et al. (2010, 2011), and Klassen et al. (2016) considered solid-body rotation. Furthermore, Meyer et al. (2018) offered a parameter scan for different initial angular velocity profiles, without turbulence. Highly supersonic turbulence leads to filamentary accretion that produces fragments with lower mass ratios to the most massive fragment, and are separated by longer distances than in a weakly turbulent cloud dominated by thermal Jeans fragmentation. The amount of turbulent fragmentation has been found to depend on the level of turbulence, with higher velocity dispersions leading to more fragmentation (Rosen et al. 2019; Fontani et al. 2018). Also, the most massive protostar moves considerably from the center of the cloud, although a (fragmenting) disk is formed in some cases (Rosen et al. 2016, 2019). The parameter scan in Meyer et al. (2018) revealed that a steep initial angular velocity profile such as the one used here produces more and earlier fragmentation.

Our choices for initial conditions are in agreement with the discussion in Beuther et al. (2018) and references therein. In that paper, observations of 20 high-mass star-forming regions showed that fragmentation is more consistent with the thermal Jeans fragmentation picture, and so low values of turbulence are expected.

Some of the studies mentioned in Sect. 9.1 consider the presence of magnetic fields. In general, the presence of magnetic fields is expected to reduce fragmentation. In Kölligan & Kuiper (2018), the effect of magnetic fields was studied in the context of high-resolution two-dimensional simulations of jet launching. Magnetic pressure was found to dominatethe inner parts of the disk, and so we expect a more stable disk in these regions. This could also mean that fragmentation may be observed later in time, as the disk would need to reach a certain size to become unstable and fragment. Magnetic braking is also expected to reduce the size of the disk, and since we observe fragmentation in the outer disk, this can reduce the number of fragments. For more details on the role of magnetic fields on forming massive stars, see for example Commerçon et al. (2011b), Seifried et al. (2011) and Rosen & Krumholz (2020).

9.4 Accretion bursts

Meyer et al. (2017, 2019) describe the existence of accretion bursts, a phenomenon we have also observed in our simulations, in the same way. Nevertheless, our analysis of the mass of the fragments (see Sect. 7.3.1), their second collapse phase, and their possible final orbits leads us to conclude that strong accretion bursts of magnitudes ≳4 as described in Meyer et al. (2019) are unlikely to happen with the predicted frequency, because the masses needed for these events correspond typically to fragments that should have undergone second collapse and, because of the reduced collisional cross section, are less probable of being accreted directly by the central massive protostar, instead forming companions.

10 Summary and conclusions

We present the highest resolution, three-dimensional self-gravity-radiation-hydrodynamical simulations of a forming and fragmenting accretion disk around a forming massive star performed to date, incorporating stellar evolution, and dust sublimation and evaporation. The high resolutions achieved in this work mean that no sink particles are necessary to describe the fragments, which also has allowed us to study the fragmentation mechanisms, interactions, and the structure of fragments. A detailed analysis of the formation, interactions, and ultimate fate of the fragments is presented. For this analysis, sophisticated data post-processing techniques allowed us to track the masses and central temperatures of the fragments and their orbits. We also studied the further evolution of these fragments into companion stars by considering the conditions for the formation of second Larson cores.

The same simulation was run five times with varying resolutions. We studied convergence in our results, and verified that we properly resolve the disk pressure scale height and the Jeans length. Our findings can be summarized as follows:

• A growing accretion disk is formed around a central massive protostar.

• The disk forms spiral arms, and they in turn fragment.

• In the fragmenting disk, the Toomre parameter indicates that the disk is globally stable while being locally unstable, and therefore an insufficiently resolved disk may appear to be Toomre-stable while fragmentation is happening at lower scales.

• Fragments are highly dynamic; they form, interact, and can be destroyed, sometimes forming secondary disks with their own spiral arms.

• Secondary disks can also fragment.

• Fragments tend to have highly eccentric and chaotic orbits, which is due to their continuous interactions. However, they canalso gain enough mass to be in more stable orbits and form companions in the middle or outer regions of the disk.

• The masses of the fragments are of the order of one solar mass once their central temperatures can reach the hydrogen dissociation limit (2000 K). Fragments that go beyond that temperature should form second Larson cores while others remain as hydrostatic cores.

• Fragments in general can migrate inwards; if they are hydrostatic cores, they might be accreted by the central massive protostar and create accretion bursts, and if they become second cores, they most likely would form spectroscopic companions.

• The number of fragments that survive the fragmentation epoch is ≲10.

• Reducing resolution by a factor of two (run x8) does not significantly affect the overall picture of the results, except for the report of lower central temperatures in the fragments.

• When reducing resolution by a factor of four (run x4), the formation of secondary spiral arms is limited, and the fragments no longer achieve the required temperatures for secondary collapse.

• When reducing the resolution by a factor of eight (run x2), the formation of the primary disk and its spiral arms is recovered, but the fragmentation physics is no longer properly achieved.

• When reducing the resolution by a factor of sixteen (run x1), only the formation and growth of the primary disk is recovered, but its internal substructure can no longer be resolved.

With each level of higher resolution another aspect of the physics of disk fragmentation appears and remains for all simulations in even higher resolution. The two highest resolution runs agree on the disk properties, its morphology, and the onset of fragmentation, while showing only subtle differences in the central temperature of the forming fragments, which suggests convergence of the physical processes involved in fragmentation. However, only the highest resolution run shows the amount of detail in the smaller scales necessary to see the formation of hydrostatic cores, that may become future companions.

Acknowledgements

G.A.O. acknowledges financial support by the Deutscher Akademischer Austauschdienst (DAAD), under the program Research Grants - Doctoral Programmes in Germany. R.K. acknowledges financial support via the Emmy Noether Research Group on Accretion Flows and Feedback in Realistic Models of Massive Star Formation funded by the German Research Foundation (DFG) under grant No. KU 2849/3-1 and KU 2849/3-2. We thank the anonymous referee for their very valuable and helpful comments and suggestions.

References

All Tables

All Figures

Fig. 1 Time evolution of the system (run x8). In the elaboration of the plots, snapshots of the simulation were taken specifically at 3.0, 8.0, and 17.5 kyr for each respective column. The rightmost column shows double the size of the other two (the dotted square corresponds to the area shown in the snapshots located to the left). For the rightmost column, notice the existence of a surviving fragment at r ~ (900, −1200) au for the midplane views.

In the text

	Fig. 2 (a) Mass, (b) luminosity, and (c) accretion rate of the central massive protostar, as formed in run x16. The dotted line in panel b corresponding to the stellar luminosity is calculated using the evolutionary tracks of Hosokawa & Omukai (2009). The gray line in panel c shows a reference accretion rate of 10−3 M⊙ yr−1.
In the text

	Fig. 3 Variation of the time-averaged mass of a sample of fragments as a function of the assumed radius taken for integration.
In the text

	Fig. 4 (a)–(e) Radial profiles of several quantities for the background disk. The intervals in profiles (a), (d), and (e) mean that the quantity has been time-averaged in that interval, excluding the endpoint. The filled diamond indicates the value at the disk radius. (f) Accretion rate into the central massive protostar as a function of time.
In the text

	Fig. 5 Comparison of the Toomre parameter values computed for the background disk (Qdisk), the mean value (Qmean), and the minimum values corresponding to the effects of the fragments and spiral arms (corresponding to Qmin). These curves were calculated with run x16, and they are time-averaged for the fragmentation epoch.
In the text

	Fig. 6 Fragment creation mechanisms. The color scale for the midplane density maps is the same as the one shown in Fig. 1.
In the text

	Fig. 7 (a) Number of fragments with a minimum lifetime of 200 yr present in the simulation as a function of time. (b) Radial position at formation time for each fragment, color-coded with the fragment lifetime.
In the text

	Fig. 8 Mass of the fragments with life span longer than 200 yr. A fixed radius of 50 au is assumed for the calculations of this plot.
In the text

	Fig. 9 Masses of the fragments with life span longer than 200 yr, for run x16. A fixed radius of 40 au is assumed for the calculations of this plot.
In the text

	Fig. 10 (a) Orbits of fragments with a life longer than 1 kyr during the time interval [6, 10] kyr, for run x8. (b) Orbits of all the fragments with a life longer than 1 kyr during all of the simulated time for run x16. The fragment identification number is located at the starting point of the orbit, and the color code is the same as that used in Figs. 8 and 9.
In the text

	Fig. 11 Fragment–fragment interactions. Panels A3 and B4: the fragment ID is shown at the starting point of the orbit. Panel A3: the time window shown in the orbits is 1.3 kyr; panel B4: this is 1.2 kyr.
In the text

	Fig. 12 Interactions of the fragments with the spiral arms, and with the forming protostar modeled as a central sink cell. The orbit in panel A4 is plotted for a window of 2 kyr, with the fragment ID marking the starting point.
In the text

Fig. 13 Secondary disk with its spiral arms around (a) fragment 11 of run x8, and (b) fragment 12 of run x16. The background map represents density in the midplane (using the same color scale as the rest of the figures), and the arrows show the comoving velocity field. The central massive protostar is located beyond the bottom of the plotted areas in both cases. The yellow circle in (b) has a radius of 3 au.

In the text

	Fig. 14 Vertical structure of fragment 12 and its close surroundings (run x16). The yellow circle has a radius of 40 au. The density uses the same color scale as the rest of the figures.
In the text

	Fig. 15 Temperature of the fragments with lifespan longer than 200 yr, for run x8. The gray box indicates the hydrogen dissociation limit.
In the text

	Fig. 16 Temperature of the fragments with lifespan longer than 200 yr, for run x16. The gray box indicates the hydrogen dissociation limit.
In the text

	Fig. 17 Fragment-destruction mechanisms.
In the text

	Fig. 18 Integrated orbit of fragment 2 of run x16 after it enters the numerical sink cell. The orange dot indicates the initial position of the fragment, and the total time integrated was 3 kyr. The sink cell (dashed circle) has a radius of 30 au.
In the text

	Fig. 19 (a) Mass and (b) luminosity of the central massive protostar, as formed in run x8.
In the text

	Fig. 20 Convergence of different quantities with resolution. The midplane density data were taken at t = 7 500 yr for all the maps in the row. The Toomre Q was time-averaged during the fragmentation epoch.
In the text

	Fig. 21 Fragment statistics for the fragment-producing runs. The charts show the fate of the fragments with a life longer than 200 yr according to the simulation, that is without taking into account the effects of the inner boundary and hydrogen dissociation, as discussed in Sect. 7.2. The estimated number of these fragments that undergo second collapse is given in the right-hand column.
In the text

	Fig. 22 Convergence of the properties of the formed protostar.
In the text

	Fig. 23 Comparison of midplane grid cell sizes in our fixed grid and other studies that use AMR. The different relevant scales of thesystem are shown as dashed lines; these were calculated as averages for the time period [6,8] kyr, and the corresponding disk radius is shown in the gray box. The number following the label for the AMR curves indicates the minimum cell size in astronomical units.
In the text