© R. Mignon-Risse et al. 2021

1 Introduction

Massive stars (>8 M_⊙) are luminous and are one of the main sources of feedback at parsec and galactic scales, especially when they explode as supernovæ (Larson 1974). The conditions in which they form remain unclear, however. This challenging problem also poses observational issues: massive stars are rare and are located far away from us (≳1 kpc) in dense regions. Theoretical difficulties also need to be overcome: in addition to magnetic fields and kinetic energy from turbulence (Tan et al. 2014), radiative energy is to be accounted for. At late stages, massive stars are also expected to drive the expansion of HII regions (Keto 2007), which adds even more complexity to the environmental conditions of their birth. The accretion mechanism in low-mass star formation is becoming increasingly understood in terms of disk-mediated accretion, but it is not yet explained for high-mass star formation. Furthermore, several outflow models rely on disk accretion, which means that a first uncertainty about the origin of outflows follows from the uncertainty attached to the accretion process. Finally, massive stars are located in stellar systems of higher multiplicity than low-mass stars (Duchêne & Kraus 2013), but the origin of this trend is unknown. It is in particular unclear whether it arises from core or disk fragmentation. These three questions about accretion, ejection, and fragmentation (as a cause for multiplicity) are tightly linked and depend on common physical ingredients: magnetic fields, turbulence, and radiation. These need to be modeled together. All of them are addressed in this suite of two papers.

Low-mass star formation is better understood than high-mass star formation, and first studies neglected the last two ingredients, namely turbulence and radiation. The first attempts to model massive star formation have therefore built on this. First of all, in the competitive accretion scenario (Bonnell et al. 2001), low- and high-mass protostars feed from a common gas reservoir after large-scale core fragmentation. In this case, the final stellar mass is not correlated to the core mass. Conversely, the scaled-up version of low-mass star formation, presented in the turbulent core accretion model (McKee & Tan 2003), proposes that massive stars form in isolation from the collapse of high-mass prestellar cores, stabilized against gravitational collapse by turbulent motions and magnetic fields. This model is hamperedby the rareness of high-mass prestellar cores (Motte et al. 2018), although candidates exist (e.g., Nony et al. 2018), and formation in isolation may not be the most common procedure. Global models such as the global hierarchical collapse model (GHC hereafter, Vázquez-Semadeni et al. 2016) and the inertial-inflow model (II hereafter, Padoan et al. 2020) challenge this core-fed accretion scenario and instead prefer large-scale dynamics, either due to collapse (GHC) or to the inertial motions following supernova feedback (II). They suggest the inclusion (and possible requirement) of turbulence. It is not clear yet whether accretion should be clump-fed or core-fed andoccur through disks or turbulent filaments (see, e.g., Rosen et al. 2019), but recent observations place (sparse but) increasingly convincing constraints on disk-mediated accretion based on unprecedented angular resolution.

Current observational constraints on disks indicate radii of 20 to 330 AU (Patel et al. 2005; Girart et al. 2018; Kraus et al. 2010) and masses of 1–8 M_⊙ (Patel et al. 2005). A complete massive star formation model should explain how these variations might be caused by the protostar environment(magnetic field strength, turbulence, or geometry). Disks can be subject to fragmentation, which might lead to the formation ofmultiple stars that are gravitationally bound. The massive hot core region G351.77-0.54 observed with ALMA at sub-40 AU resolution reveals 12 substructures within a few thousand AU with a broad range of core separations (Beuther et al. 2019), consistent with thermal Jeans fragmentation of a dense core, and possibly with the global hierarchical model (Vázquez-Semadeni et al. 2016). The disk in HH 80–81 might be prone to fragmentation (Fernández-López et al. 2011)as well. There does not seem to be any general trend of the disk stability, but the advent of ALMA will increasethe statistics. In the high-mass star-forming region IRAS 23033+5951, four millimeter sources are identified with the Northern Extended Millimeter Array (NOEMA) and the IRAM telescope, two of which exhibit protostellar activity. One of them is stable and the other is prone to fragmentation in the inner 2000 AU (Bosco et al. 2019). Disk fragmentation may either lead to the formation of companion stars or to the accretion of clumps onto the central object. Moreover, accretion leads to a radiative shock at the stellar surface and the energy is radiated away, therefore clump accretion can be detected by luminosity outbursts. The process of this accretion luminosity is not fully understood, in particular, the conditions under which the radiation would escape rather than being advected together with the gas. It appears to be a recurrent mechanism in low-mass star formation, however, and relies on disk-mediated accretion and star–disk interaction. It therefore also advocates the same accretion method for the formation of high-mass stars (Caratti o Garatti et al. 2017).

Despite the lack of systematic constraints, the presence of disks around young massive protostars (L < 10⁵ L_⊙) is well established in general (see the reviews by Beltrán & de Wit 2016; Beltrán 2020). Disk properties may set the initial conditions for the formation of multiple stellar systems, and they strongly depend on the threading magnetic fields.

Constraints on magnetic field structures and strength are recent and have been acquired with new polarimetric instruments. In a sample of 21 high-mass star-forming clumps, subparsec magnetic fields appear to be structured (Zhang et al. 2014). The hourglass shape that arises because the field lines are pulled by the collapsing gas is present (Beltrán et al. 2019), as in low-mass protostellar systems (e.g. Maury et al. 2018). The parameter $\mu =(M/\varphi )/{(M/\varphi )}_{\text{crit}}$ is the ratio of mass-to-flux to critical mass-to-flux, where ϕ is the magnetic flux. It indicates whether magnetic fields can (μ < 1) or cannot (μ > 1) prevent collapse on their own. Nevertheless, a μ ≳ 1 still affects the gas dynamics. Several studies agree on supercritical values of μ = 1–4 (Falgarone et al. 2008) or even μ ~ 1–2 (Girart et al. 2009; Li et al. 2015; Pillai et al. 2015), suggesting an important role of magnetic fields. Quantitatively, the field strength has the order of 0.1–1 mG in a sample of infrared dark clouds (IRDCs, Pillai et al. 2016) and in an ultracompact HII region (UCHII, Tang et al. 2009), based on the Chandrasekhar-Fermi method. The Chandrasekhar–Fermi method relates the plane-of-the-sky field strength with the line-of-sight velocity dispersion, using the phase velocity of transverse Alfvén waves. (Chandrasekhar & Fermi 1953). Comparisons with magnetohydrodynamical (MHD) simulations have shown that fragmentation is consistent with turbulence dominating the magnetic energy (Palau et al. 2013; Fontani et al. 2016). Nonetheless, magnetic energy has been found to be comparable to (Falgarone et al. 2008; Girart et al. 2013) or to dominate the turbulent energy (sub-Alfvénic turbulence, Pillai et al. 2015) in several sources. Girart et al. (2013) have found equipartition of rotational energy, magnetic energy, and turbulent energy in a fast-rotating core, with μ = 6, indicating three mechanisms capable of slowing down the collapse.

This shows that observations agree on the presence of disk-like structures, with several occurrences of fragmentation and sizes of tens to several hundred AU. Magnetic fields have non-negligible strengths, and their presence may affect disk formation (and subsequently, outflow launching) and its fragmentation into multiple stellar systems. We summarize the recent improvements in our understanding of disk formation on the side of numerical studies, focusing on the physics they have included, and in particular, the treatment of MHD and radiative transfer.

Disk-mediated accretion for massive protostars has emerged in multidimensional simulations as part of the so-called flashlight effect (Yorke & Sonnhalter 2002; Kuiper et al. 2010a) to overcome the radiation barrier problem (Larson & Starrfield 1971). This effect describes the disk thermal radiation as being radiated preferentially off the plane, ending in a small radiative force against the accretion flow. This is different from the unidimensional view. Progress has also been made in the low-mass star formation context with the inclusion of magnetic fields in numerical simulations in the ideal MHD frame (e.g., Fromang et al. 2006). Many studies have shown that the flux-freezing condition in a collapsing core leads to the accumulation of magnetic fields in the central region, which induces a strong magnetic braking and prevents disk formation (see, e.g. Hennebelle & Fromang 2008, and Seifried et al. 2011 in the high-mass regime). This is referred to as the magnetic catastrophe because many disks are observed around low- and high-mass protostars (Cesaroni et al. 2005). Three ingredients have been introduced and shown separately to allow the formation of disks and reconcile numerical simulations and observations in this respect: misalignment of the rotation axis and the magnetic field axis, turbulence, and nonideal MHD effects. Misalignment (Joos et al. 2012) and turbulence (Joos et al. 2013; Lam et al. 2019 for low-mass stars, Seifried et al. 2012 for high-mass stars) directly reduce the magnetic braking efficiency. Nonideal (also called resistive) MHD effects, namely ambipolar diffusion (AD), Ohmic dissipation, and the Hall effect, provide a mechanism to limit the accumulation of the magnetic field strength and therefore magnetic braking. AD is probably the most frequently studied nonideal MHD effect because it starts to dominate at lower densities than the others, and indeed promotes disk formation in the low- (Masson et al. 2016)and high-mass regime (Commerçon et al. 2021, hereafter C21). In several studies, nonideal MHD appears to be the main regulator of disk formation (AD in Hennebelle et al. 2016a) even when subsonic turbulence is included (Wurster & Lewis 2020).

In parallel, most numerical studies of massive star formation have focused on the radiative transfer aspect because of the radiation pressure barrier (Larson & Starrfield 1971) and neglected magnetic fields. First radiation-hydrodynamical implementations have relied on the flux-limited diffusion (FLD) approximation (Levermore & Pomraning 1981) to describe infrared radiation interacting with dust-gas mixture. The FLD method is well suited for radiation transport in optically thick media, but is not adapted to strongly anisotropic radiation fields. The particular treatment of stellar radiation, also called irradiation, was later improved to track the higher-energy stellar photons and accurately compute the opacity during the first interaction of the photons with the ambient medium (Kuiper et al. 2010b; Flock et al. 2013; Ramsey & Dullemond 2015; Rosen et al. 2017; Mignon-Risse et al. 2020; Gressel et al. 2020; Fuksman et al. 2021).

The common inclusion of radiative transfer and MHD in numerical codes has shown that both effects contribute to limiting fragmentation. Commerçon et al. (2011a) showed the prevention of early core fragmentation, while Myers et al. (2013) obtained similar results at later times. Secondary fragmentation, responsible for the formation of companion stars, is also inhibited, as found by Peters et al. (2011).

As presented above, the modeling of magnetized disks requires nonideal MHD effects to circumvent the so-called magnetic catastrophe (C21). In this work, we extend the study by C21 and present the first numerical simulations that include a hybrid radiative transfer method and nonideal MHD (ambipolar diffusion) with the aim to identify the accretion conditions of massive protostars, and their outflow-launching mechanism (Mignon-Risse et al, in prep., hereafter Paper II) with realistic physical ingredients. To do so, we consider an initial velocity field consistent with turbulence (of various amplitudes, corresponding to several runs) in order to mimic nonidealized environmental conditions for the birth of a massive protostar.

This study is organized as follows. The numerical methods are presented in Sect. 2. In Sect. 3, we analyze the disk-mediated accretion, emphasizing the primary, secondary, and circumbinary disks properties and the alignment of disk and magnetic field. Our results and the limitations of the study are discussed in Sect. 4. We conclude in Sect. 5.

2 Methods

In this section, we present the set of equations that is solved numerically, and the set of initial conditions. We summarize our sink particle algorithm, and finally, we present the physical criteria that define a disk in the post-processing step. This work is intended to extend the study of C21, which was made with the flux-limited diffusion method for radiative transfer, but this time, with a hybrid radiative transfer method and in a turbulent medium. An additional difference lies in the sublimation model we take here, and in the optically thin sink volume (see below).

2.1 Radiation magnetohydrodynamical model

We integrated the equations of radiation MHD in the adaptive-mesh refinement (AMR) RAMSES code (Teyssier 2002; Fromang et al. 2006) with ambipolar diffusion (Masson et al. 2012), and the so-called hybrid radiative transfer method (Mignon-Risse et al. 2020), namely the M1 method (Levermore 1984; Rosdahl et al. 2013; Rosdahl & Teyssier 2015) for stellar radiation and the flux-limited diffusion (FLD, Levermore & Pomraning 1981; Commerçon et al. 2011b; Commerçon et al. 2014)otherwise. The set of equations we solved is $$\begin{array}{cc}\frac{\partial \rho }{\partial t}+\nabla \cdot [\rho u]=& 0,\\ \frac{\partial \rho u}{\partial t}+\nabla \cdot [\rho u\otimes u+P\mathbb{I}]=& -\lambda \nabla E_{\text{fld}}+\frac{{\kappa }_{\text{P},\star }\rho }{\text{c}}{F}_{\text{M1}}+F_{\text{L}}-\rho \nabla \varphi ,\\ \frac{\partial E_{\text{T}}}{\partial t}+\nabla \cdot [u\left(E_{\text{T}}+P+B^2/2\right)& \\ -(u\cdot B)B-E_{\text{AD}}\times B]=& -{\mathbb{P}}_{\text{fld}}\nabla :u+{\kappa }_{\text{P},\star }\rho \text{c}{E}_{\text{M1}}-\lambda u\nabla E_{\text{r}}\\ & +\nabla \cdot \left(\frac{\text{c}\lambda }{\rho {\kappa }_{\text{R},\text{fld}}}\nabla E_{\text{r}}\right)-\rho u\cdot \nabla \varphi ,\\ \frac{\partial E_{\text{M1}}}{\partial t}+\nabla \cdot F_{\text{M1}}=& -{\kappa }_{\text{P},\star }\rho \text{c}{E}_{\text{M1}}+{\dot{E}}_{\text{M1}}^{\star },\\ \frac{\partial F_{\text{M1}}}{\partial t}+{\text{c}}^2\nabla \cdot {\mathbb{P}}_{\text{M1}}=& -{\kappa }_{\text{P},\star }\rho \text{c}{F}_{\text{M1}},\\ \frac{\partial E_{\text{fld}}}{\partial t}-\nabla \cdot \left(\frac{\text{c}\lambda }{\rho {\kappa }_{\text{R},\text{fld}}}\nabla E_{\text{fld}}\right)=& {\kappa }_{\text{P},\text{fld}}\rho \text{c}\left({\text{a}}_{\text{R}}{T}^4-E_{\text{fld}}\right),\\ \frac{\partial B}{\partial t}-\nabla \times \left[u\times B+E_{\text{AD}}\right]=& 0,\\ \nabla \cdot B=& 0,\\ \Delta \varphi =& 4\pi \text{G}\rho .\end{array}$$(1)

Here, ρ is the material density, u is the velocity, P is the thermal pressure, λ is the FLD flux limiter, E_fld is the FLD radiative energy, κ_P,⋆ is the Planck mean opacity at the stellar temperature, F_M1 is the M1 radiative flux, F_L = (∇ × B) × B is the Lorentz force, ϕ is the gravitational potential, E_T is the total energy E_T = ρϵ + 1∕2ρu² + 1∕2B² + E_fld (ϵ is the specific internal energy), E_M1 is the M1 radiative energy, B is the magnetic field, E_AD is the ambipolar electromotive force, ${\mathbb{P}}_{\text{fld}}$ is the FLD radiative pressure, κ_P,fld is the Planck mean opacity in the FLD module, κ_R,fld is the Rosseland mean opacity, a_R is the radiation constant, ${\mathbb{P}}_{\text{M1}}$ is the M1 radiative pressure, and ${\dot{E}}_{\text{M1}}^{\star }$ is the stellar radiation injection term.

We injected radiative energy into the M1 module only with the primary sink, while other sinks (formed after the primary one) radiated within the FLD module. The justification for this is twofold. First, in the M1 method, the radiative fluxes from several sources are added, whereas two radiation beams should not interact when they cross each other (González et al. 2007). The generalization of the hybrid method, which relies on the M1 method, to several stellar sources should be addressed in dedicated studies. Second, we address the origin of outflows around individual massive protostars (paper II), and using the M1 method for treating the radiation of one star is sufficient for this. The ambipolar electromotive force is equal to $$E_{\text{AD}}=\frac{{\eta }_{\text{AD}}}{B^2}\left[(\nabla \times B)\times B\right]\times B,$$(2)

where η_AD is the ambipolar diffusion resistivity. It depends on the density, temperature, and magnetic field strength. The resistivities were precomputed using a chemical network to calculate the equilibrium abundances of the molecules (neutrals) and main charge carriers in conditions of prestellar core collapse (Marchand et al. 2016), depending on the density and temperature. Chemical equilibrium was assumed because the associated timescale is shorter than the free-fall time at the densities we considered.

The term κ_P,⋆ρcE_M1 couples the M1 and the FLD methods through the equation of the evolution of the internal energy, $$C_{\text{v}}\frac{\partial T}{\partial t}={\kappa }_{\text{P},\star }\rho \text{c}{E}_{\text{M1}}+{\kappa }_{\text{P},\text{fld}}\rho \text{c}\left(E_{\text{fld}}-{\text{a}}_{\text{R}}{T}^4\right).$$(3)

We used the ideal gas relation for the internal specific energy ρϵ = C_vT, where C_v is the heat capacity at constant volume.

2.2 Physical setup

We assumed similar initial conditions as C21. The free-fall time in the central plateau (which contains ~ 15 M_⊙) of the density profile is then $${\tau }_{\text{ff}}=\sqrt{\frac{3\pi }{32\text{G Bar}{\rho }_{\text{0}}}}\simeq 24\text{kyr},$$(4)

where G is the gravitational constant, and $\text{Bar}{\rho }_{\text{0}}$ is the density of the central plateau.

The initial core was threaded by a uniform magnetic field oriented along the x-axis. We set the magnetic field strength by the ratio of mass-to-flux to critical mass-to-flux, $\mu ={(M/\Phi )}_{\text{0}}/{(M/\Phi )}_{\text{crit}}$, where ${(M/\Phi )}_{\text{0}}=M_{\text{c}}/(\pi B_{\text{0}}{R}_{\text{c}}^2)$ and ${(M/\Phi )}_{\text{crit}}=0.53/(3\pi )\sqrt{5/G}$ (Mouschovias & Spitzer 1976). Strong (μ = 2, B₀ = 170 μG) and moderate (μ = 5, B₀ = 68 μG) magnetic fields were considered. A drawback of this uniform distribution is that the mass-to-flux ratio decreases as μ ~ 1∕R and is higher in the inner parts of the core, with μ ≈ 50 in the central plateau (for runs with μ = 5). This corresponds to a weakly magnetized medium. We expect the central magnetic field strength to increase as B ∝ ρ^2∕3 as the core contracts and before ambipolar diffusion starts to dominate (at ρ ~ 10⁻¹⁵g cm⁻³), however. The magnetic field therefore plays a dynamical role in the collapse (see C21).

An initial velocity dispersion was imposed to mimic a turbulent medium, and it followed a Kolmogorov power spectrum P(k) ∝ k^−5∕3, similar to Commerçon et al. (2011a). Phases were randomly sampled, and one realization was considered, that is, we did not vary the velocity field, but only its amplitude between two runs. The turbulence was not sustained, but the turbulent crossing time (~ 0.5–2 Myr with T = 20 K, see below) was significantly longer than the simulation time here, so that it is not expected to affect our results. A low level of (solid-body) rotation, E_rot ∕E_grav = 1%, was initially imposed around the x-axis and dominated the specific angular momentum in subsonic runs. We considered four runs (see Table 1) in which we varied the initial Mach number $\mathcal{M}$ and Alfvénic Mach number ${\mathcal{M}}_{\text{A}}$. The runs can be divided into two classes that we refer to throughout the paper, depending on the relative effect of turbulence and magnetic fields. Turbulence is sub-Alfvénic in runs NOTURB (${\mathcal{M}}_{\text{A}}=0$) and SUBA (${\mathcal{M}}_{\text{A}}<1$), and super-Alfvénic in runs SUPA and SUPAS (${\mathcal{M}}_{\text{A}}>1$). Regarding the Mach number, runs SUPA and SUBA have subsonic turbulence with $\mathcal{M}=0.5,$ while run SUPAS has supersonic turbulence with $\mathcal{M}=2$.

Sink particles are introduced to mimic the presence of a protostar and radiate as blackbodies. Their radii and internal luminosities (which give their effective surface temperature) are interpolated from a pre-main-sequence track (Kuiper & Yorke 2013), based on their mass at a given time and their accretion rate averaged over their lifetime. We modeled the dust sublimation by decreasing the dust-to-gas ratio progressively with the temperature. Gray opacities were taken from Semenov et al. (2003), and the dust sublimation was modeled with a dust-to-gas ratio that vanished at high temperatures, similarly to Kuiper et al. (2010a). The dust-to-gas ratio varies as $$\frac{M_{\text{dust}}}{M_{\text{gas}}}(\rho ,T)={\left(\frac{M_{\text{dust}}}{M_{\text{gas}}}\right)}_0\left(0.5-\frac{1}{\pi }\text{arctan}\left(\frac{T-T_{\text{evap}}(\rho )}{100}\right)\right),$$(5)

where ${\left(\frac{M_{\text{dust}}}{M_{\text{gas}}}\right)}_0=1\%$ is the initial dust-to-gas mass ratio, and the evaporation temperature is $$T_{\text{evap}}(\rho )=g{\left(\frac{\rho }{1\text{g}{\text{cm}}^{\text{-3}}}\right)}^{\beta },$$(6)

where g = 2000 K and β = 0.0195 (Isella & Natta 2005). The gas opacity was taken to be equal to 0.01 cm² g⁻¹, so that the total opacity tends toward this value as the temperature increases beyond T_evap. Finally, we set the opacity in the primary sink particle volume to a value chosen so that the local optical depth was the minimal optical depth allowed numerically (10⁻⁴). This floor value for the optical depth is a numerical parameter used in optically thin cells in order to gain performance with the FLD solver (see Appendix A of Vaytet et al. 2018). Gas and radiation are hardly modeled within the sink volume, but stellar radiation is meant to escape this volume. This justifies our subgrid model for decoupling gas and radiation within the sink volume (see the discussion in Appendix A).

2.3 Resolution and sink particles

Boundary conditions are periodic and the simulation box is 0.8 pc large, hence the gravitational effects due to the periodicity are marginal. Indeed, the box length is twice the core diameter. At the core border, the gravitational acceleration exerted on a gas particle scales as $a_{\text{grav}}~M_{\text{c}}/R_{\text{c}}^2$. In comparison, the acceleration due to the nearest (0.6 pc) periodic core is $~M_{\text{c}}/{(3R_{\text{c}})}^2=a_{\text{grav}}/9$. A cell effective resolution is 0.5^ℓ (in units of the box length), where ℓ is the AMR level of refinement. The coarse grid is 32³ (level ℓ = 5), with 10 additional levels of refinement. It leads to a smallest cell width of 5 AU. We run a similar set of simulations with a maximal resolution of 10 AU to probe resolution convergence. We use the prefix LR, as “low-resolution”, to refer to these runs (not shown here for conciseness). Refinement is performed on the Jeans length: it must be sampled by at least 12 cells, following Truelove et al. (1997). Sinks were introduced at the finest AMR level, after a clump has been found to be bound and Jeans unstable (see C21, Bleuler & Teyssier 2014). Sinks accrete the environmental material that enters their accretion volume as follows. Theradius of this sphere, the so-called accretion radius, was set to four times the finest resolution, that is, 20 AU (40 AU for the LR runs). We used a threshold accretion scheme based on the Jeans density in order to avoid artificial fragmentation. If the density of a sink cell exceeds the local Jeans density, 10% of the excess is accreted by the sink particle within one time step. All of the accreted mass is added to the sink particle mass (no outflow subgrid model). Sink particles are collisionless, they interact only gravitationally with the gas and with their companions. Their gravitational potential was prevented from diverging by the use of a so-called softening length, which was set equal to the accretion radius. Finally, sink particles could merge if their accretion radii overlapped. We did not add any criterion to prevent sinks from merging.

2.4 Disk identification

The primary disk properties are presented in Sect. 3.3. They were computed from the cell-by-cell disk selection, which relies on the following criteria (presented in Joos et al. 2012). It is rotationally supported: $\rho v_{\varphi }^2/2>f_{\text{thres}}P$, where v_ϕ is the azimuthal velocity and P is the thermal pressure. We chose f_thres = 2 as in Joos et al. (2012). The gas number density is higher than n = 10⁹cm⁻³, i.e. ρ = 3.85 × 10⁻¹⁵g cm⁻³. The gas is not about to free-fall onto the central object too rapidly: v_ϕ > f_thresv_r, where v_r is the radial velocity. The disk vertical structure is in hydrostatic equilibrium: v_ϕ > f_thresv_z, where v_z is the vertical velocity. We note that there is no connectivity criterion, therefore the extremity of a high-density spiral arm can be considered to be part of the disk while the inter-arm low-density region (due to the gas being swept) may not be. We computed the disk radius as the radius enclosing 90% of the total disk selection mass to avoid accounting for transient negligible contributions from larger scales. We added a geometrical criterion to avoid perturbations from companions and their possible disks: the cell must be located at less than 0.9 times the binary separation. This is justified a posteriori as individual disks are observed around each star.

3 Results

We first summarize some of the main and common features of all runs in Sect. 3.1. The sink mass evolution, which depends on disk-mediated accretion, is presented in Sect. 3.2. Then we focus on the properties of the primary disk (Sect. 3.3), the secondary disk (Sect. 3.4), and the circumbinary disk (Sect. 3.5). Finally, we study the primary disk alignment with the core-scale magnetic fields (Sect. 3.6).

3.1 Overview and common features

3.1.1 Formation of structures and stars

In the four simulations, as the gravitationally unstable cloud core collapses, the first sink particles form at t ≈ 29 kyr. Figure 1 shows density slices in the disk plane and perpendicular to the disk plane in the four runs at time t = 50 kyr, together with gas velocity and magnetic field lines. Except in run SUPAS (the most turbulent, third column of Fig. 1), where a large filament-like structure forms due to the stronger inner turbulent support (see below), the dense region (ρ ≳ 10⁻¹⁶g cm⁻³) rapidly concentrates in a sphere of diameter ~2000 AU (center panels, second row of Fig. 1). This is reminiscent of the structure described by Machida & Hosokawa (2020), and we attributed itto the inability of the toroidal magnetic pressure to launch an outflow because of turbulence. We show in Fig. 2 the plasma beta (β = P_th∕P_mag, where P_th and P_mag are the thermal pressure and magnetic pressure, respectively) corresponding to Fig. 1. The central region in run SUPAS (third column, second row) is magnetically dominated (β < 1) at 2000 AU scales while thermally dominated (β > 1) matter is infalling. This illustrates the importance of accurately accounting for the coupling of gas and magnetic fields in order to assess the dynamical role expected to be played by magnetic effects. Accretion disks form in all runs around the primary sink (third and last rows of Fig. 1). In runs NOTURB and SUBA, in which turbulence is sub-Alfvénic, no secondary sink forms. With super-Alfvénic turbulence (runs SUPA and SUPAS), a secondary long-lived sink particleforms in the primary sink accretion disk. We study the stellar multiplicity in Sect. 3.2 and the secondary and circumbinary disk properties in Sect. 3.4 and Sect. 3.5 in more detail. When not mentioned otherwise, we refer to the primary sink and to the primary disk for conciseness.

3.1.2 Magnetic field evolution

As collapse occurs, the magnetic field strength is expected to increase in the central regions. The histograms of the density-magnetic field strength for the different runs are shown in Fig. 3. At densities below ~ 10⁻¹⁵ g cm⁻³, we recover the ideal MHD limit within which B increases with ρ. In runs NOTURB and SUBA, the high-B, low-ρ part of the histogram is populated by outflowing material ejected from the most magnetized regions. At high densities, the plateau-like feature is present in the four runs and contrasts with ideal MHD calculations, as shown in the low-mass (Masson et al. 2016) and high-mass (C21) regimes. This is due to ambipolar diffusion, which becomes dominant above ρ ≳ 10⁻¹⁵ g cm⁻³. The diffusion coefficient varies nonlinearly with the magnetic field strength η_AD ∝ B²∕ρ, which explains its strong regulating effect. The plateau is located between ~ 0.1 G in the super-Alfvénic runs (SUPA, SUPAS) and at ~0.3 G in the sub-Alfvénic runs (NOTURB, SUBA). The inclusion of ambipolar diffusion prevents the magnetic field strength from increasing, which would change the disk structure and possibly the outflows because a strong magnetic field is reasonably expected in the magnetocentrifugal mechanism (C21). The large dispersion observed in the bottom left panel of Fig. 3 at low density is caused by turbulence in the core.

3.1.3 Asymmetries

While the numerical setup in run NOTURB is initially axisymmetric and symmetric with respect to the (x = 0, y, z)-plane (that we refer to as the “north-south" symmetry), these symmetries are broken. First, pockets of magnetized plasma (β < 1) are regularly expelled from the disk outer edge (top panels of Fig. B.1). This is visible in run NOTURB, but hardly seen in the other turbulent runs. We investigate in Appendix B whether the magnetic interchange instability, which has been found to redistribute magnetic flux after accumulation around sink particles (Krasnopolsky et al. 2012) or at the ambipolar diffusion radius (Li et al. 2011), is responsible for this. The timescale associated with the interchange instability is indeed short enough (compared to the local advection timescale) for the interchange instability to be a good candidate. This phenomenon is not the only asymmetry arising in the simulation. Filamentarystructures link the densest regions (where the sink-disk system is) to the envelope, which we refer to as “streamers”. They are visible in Fig. 2, as the filaments that are dominated by thermal pressure (β > 1), unlike the gas that surrounds them. They have a density ρ ≳ 10⁻¹⁵g cm⁻³. They appear to be a path for the accretion flow and pull the magnetic field lines along with them, which in turn form an hourglass shape. Because the Lorentz force has no component parallel to the field lines, the gas can move along the lines to join the streamers without any magnetic resistance, in a similar way as the bead-on-a-wire picture for magnetocentrifugal jets. These streamers form perpendicular to the core-scale magnetic field in all runs. In run NOTURB, this plane is also that of the accretion disk. Nonetheless, they connect to the disk outside of the disk plane (first column, third and last rows of Fig. 2) either from above or below, breaking the north-south symmetry. We attribute this to the strong magnetic forces around the streamers. In runs SUPA and SUBA, the streamers are much thicker. This might be an indication of turbulent diffusion, but is not investigated further here. This gives rise to the filament-like structure of width ~ 2000 AU in run SUPAS, as mentioned above. Overall, the streamers do not appear to set the disk formation plane (studied in Sect. 3.6). Nevertheless, the symmetry breaking they provide is important to us, considering that 16% of the outflowsreported in Wu et al. (2004) are monopolar, and our aim is to study a nonidealized case that could be compared to observations.

Fig. 1 Density slices perpendicular (first and second row, zoomed by a factor of 10) and parallel (third and fourth row, zoomed by a factor of 5) to the disk plane, at t ≈ 50 kyr. Streamlines corresponding to magnetic field lines and arrows corresponding to the velocity field are overplotted. Columns from left to right: Runs NOTURB, SUPA, SUPAS, and SUBA. The mass density ρ = 10−19 g cm−3 corresponds to the particle density n = 2.6 × 104 cm−3 and ρ = 10−11 g cm−3 to n = 2.6 × 1012 cm−3. White dots represent sink particles.

Fig. 2 Same as Fig. 1, but with plasma beta (ratio of the thermal pressure and the magnetic pressure) slices. The gas moves along magnetic field lines until it forms thermally dominated infalling filaments. Disks have β > 1 as well.

3.1.4 Angular momentum transport

We report three mechanisms for angular momentum transport to help accretion onto the central object. First, outflows are ubiquitous except in run SUPAS, where a monopolar and transient outflow develops (see Paper II). Second, spiral arms are present in the disk (similar to Klassen et al. 2016) and exert gravitational torques that transport angular momentum outward. Third, magnetic braking occurs. We computed contributions to the angular momentum flux as (Joos et al. 2012) $$F_{\text{out}}=|{\displaystyle {\int }_{\text{S}}\rho }r{v}_{\varphi }v\cdot \text{d}S|$$(7)

for the outflows (using the selection criteria presented in Paper II), $$F_{\text{grav}}^{\text{r}}=|{\displaystyle {\int }_{S_{\text{rad}}}r}\frac{g_{\varphi }}{4\pi \text{G}}g\cdot \text{d}S|$$(8)

for the gravitational torque, and $$F_{\text{mag}}=|{\displaystyle {\int }_{\text{S}}r}\frac{B_{\varphi }}{4\pi }B\cdot \text{d}S|$$(9)

for the magnetic torque. The first and third integrals were performed over the surface S of a cylinder centered onto the primary sink, of radius R = 1000 AU and height H = 1000 AU, and they were oriented along the angular momentum vector. The surface S_rad only accounts for the surfaces of this cylinder that are perpendicular to the cylindrical radial vector because we expect gravitational torques to transport angular momentum radially rather than vertically. In order to have comparable values from one run and one time step to the next, we divided these fluxes by the mass within the cylinder. Figure 4 shows F_out∕M, F_mag∕M, and $F_{\text{grav}}^{\text{r}}/M$ as a function of the sink age. Angular momentum transported by outflows is slightly higher in run NOTURB than in SUPA and SUBA. We find higher angular momentum transport from magnetic torques than from outflows. Magnetic braking is initially stronger in the initially most magnetized model (run SUBA). After a sink age ~ 20 kyr, it is lower inrun SUBA than in NOTURB, suggesting that the turbulence reduces magnetic braking. This is confirmed by the even lower magnetic braking in runs SUPA and SUPAS. In the right panel of Fig. 4, the gravitational torque is stronger with increasing turbulence. The magnetic torque generally dominates the angular momentum transport, except in run SUPAS at later times, where magnetic and gravitational torques are comparable. We performed the same comparison with R = 100 AU, and the results remain qualitatively unchanged.

Fig. 3 Density-magnetic field strength histograms at t = 50 kyr. From left to right: run NOTURB and SUPA (top), SUPAS and SUBA (bottom).

Fig. 4 Evolution of angular momentum transported by the outflows (left panel), magnetic fields (middle panel), and by gravitation (right panel) within a cylinder of radius and height 1000 AU.

3.2 Sink mass history and stellar multiplicity

We first studied the mass evolution of the sink particles. The left panel of Fig. 5 displays the primary and secondary (when there is one) sink mass as a function of primary sink age. The most massive star formed is M ≃ 15.8 M_⊙ in run NOTURB, at the end of the run, when the primary sink age is 55 kyr. Globally, two different behaviors are visible in the sub-Alfvénic (NOTURB, SUBA) and super-Alfvénic runs (SUPA, SUPAS). There is a delay of ~8 kyr between runs SUPA and SUPAS, but a comparable slope (mean accretion rate). The mass accretion is much smoother in the sub-Alfvénic cases than in the super-Alfvénic cases. This is confirmed by the accretion rate, displayed in the right panel of Fig. 5 for runs NOTURB and SUPA. It is smoothed over a temporal period of ~ 1 kyr. The values for runs SUBA and SUPAS are not displayed here for readability and show similar features to runs NOTURB and SUPA. The values are mainly between 10⁻⁴ and 10⁻³ M_⊙ yr⁻¹ in run NOTURB. The accretion rate over the primary sink in run SUPA, which includes initial turbulence, is first comparable to NOTURB. After ~12 kyr, it becomes erratic. Instantaneously (i.e., without the smoothing), it has zero values most of the time. We recall that our sink accretion scheme relies on a high-enough density and Jeans-unstable gas within the sink volume. The absence of accretion at a given time therefore means that the sink volume has not gathered enough material to be accreted. The main accretion events in the super-Alfvénic runs are more dramatic than in the sub-Alfvénic runs, with companion sink particles or orbiting massive clumps increasing the primary sink mass by a fraction of a solar mass instantaneously. Averaged over time, we obtain accretion rates that agree with observational values (Motte et al. 2018 and references therein).

We report the formation of a long-lived binary system in the two super-Alfvénic runs (SUPA and SUPAS, see Table 2). In run SUPAS, three additional sink particles form from initial fragmentation and four form in the disk plane, but merge with the primary or secondary sinks. The secondary sink forms ≈ 17 kyr after the primary.This occurs at the extremity of a spiral arm.The secondary particle survives until the end of the run, that is, after a lifetime ≳ 33 kyr. This systemalso forms in the corresponding lower-resolution run (LRSUPAS), at an age difference of ≈ 18 kyr instead of ≈ 17 kyr. Because the resolution was lower, run LRSUPAS was carried out up to t ~ 106 kyr, and the secondary sink particleis ≈60 kyr old. The stellar masses are then 9.8 M_⊙ and 8.9 M_⊙. The same formation mechanism leads to the birth of a long-lived companion in run SUPA at a primary sink age of ≈ 19 kyr. Similarly, a binary system forms in LRSUPA as well, but at later times (≈ 43 kyr after the primary), with final masses of 19.7 M_⊙ and 8 M_⊙ at t ≈ 102 kyr. The secondary sink is 30 kyr old. Interestingly, after t ~ 78 kyr, the primary sink gains only a fraction of a solar mass within more than 20 kyr, while the companion accretes 6 M_⊙. Overall, we obtain long-lived (at least tens of kyr) binary systems in the super-Alfvénic runs. This demonstrates the effect of turbulence on fragmentation even when magnetic fields are present at a moderate level.

The sinks forming the binary system in runs SUPA and SUPAS have mass ratios of about unity (in run LRSUPA where it may tend toward 1) and even greater than1 in run SUPAS. The sink mass history shows (left panel of Fig. 5) that the primary sink is partially starved due to the presenceof the secondary sink, as compared to the sub-Alfvénic runs, in which no binary forms. In both runs, the secondary sink mass quickly becomes comparable to (and even greater than, in run SUPAS) the primary sink mass. The evolution of the secondary sink mass is very similar in runs SUPA and SUPAS, and the accretion rate is greater than for the primary sink (right panel of Fig. 5). This suggests that the primary sink accretion disk might be a more favorable place for accretion than the central location of the primary sink (when it dominates the total sink mass). The radial gravito-centrifugal equilibrium in the disk likely reduces the accretion rate onto the primary sink, but not onto the secondary sink. As mentioned above, the binary systems form from disk fragmentation (see Sect. 3.3.4) rather than core fragmentation. All the sink particles that formed from initial core fragmentation have merged. We recall that we used AMR based on the thermal Jeans length. The numerical convergence we performed with LR runs supports a physical rather than numerical fragmentation. Nevertheless, the absence of a criterion to prevent sink merging in our simulations means that the final sink number may be higher. Our multiplicity results can therefore be considered as lower-limit values to first order.

Figure 6 shows the relation of the primary sink and the disk mass (left panel; the proper disk mass evolution is displayed in the right panel and is discussed in Sect. 3.3.1). It appears that there is a correlation between several sink mass gain events and disk mass loss events in the four runs, but the masses involved are much lower in the sub-Alfvénic runs. This is consistent with the gas falling smoothly onto the central star through the accretion disk in sub-Alfvénic runs, while in super-Alfvénic runs, clumps form in the disk and are subsequently accreted.

Fig. 5 Primary and secondary sink mass (left panel) and accretion rate (right panel) as a function of primary sink age for the four runs. The accretion rate is plotted for one sub-Alfvénic (NOTURB)and one super-Alfvénic (SUPA) run for readability. It has been smoothed over periods of ~ 1 kyr.

3.3 Primary disk properties

3.3.1 Mass and radius

The disk mass temporal evolution is displayed in the right panel of Fig. 6. It globally increases with sink age, with more variations in the super-Alfvénic runs. We obtain disk masses ranging from ≈ 1–8 M_⊙ (for t > 10 kyr). This confirms the trend observed in C21 and extends it to a turbulent medium: in the hydrodynamical case, disks are more massive (10 M_⊙) than in the presence of magnetic fields.

Figure 7 shows the primary disk radius as a function of time (left panel) and its comparison with the analytical prediction of Hennebelle et al. (2016a). As shown in the left panel of Fig. 7, the disk radius is between 50 and 200 AU in all runs most of the time. Large and pointlike increases coincide with the presence of a large spiral arm. The quasiperiodic variations found in runs SUPA and SUPAS are due to the orbital motions that affect the gas dynamics between the two stars.

Fig. 6 Primary sink mass against the disk mass (left panel) and disk mass as a function of time (right panel) for the four runs. In the right panel, colored circles indicate the secondary sink formation epoch.

Fig. 7 Disk radius (left panel) and ratio of the disk radius and the theoretical value (right panel, Eq. (10)) as a function of time for the four runs. Colored circles indicate the secondary sink formation epoch.

3.3.2 Semianalytical estimate of the disk radius

We compared the disk sizes with the theoretical predictions from Hennebelle et al. (2016a) for magnetically regulated disk formation with ambipolar diffusion. They were obtained from the equality between various timescales at the centrifugal radius: On the one hand, the timescale required to generate a toroidal field from differential rotation and the timescale required for ambipolar diffusion to diffuse it vertically. On the other hand, magnetic braking and rotation timescales were set equal as well. The disk radius set by ambipolar diffusion is then $$r_{\text{d},\text{AD}}\simeq 18\text{AU}\times {\delta }^{2/9}{\left(\frac{{\eta }_{\text{AD}}}{0.1\text{s}}\right)}^{2/9}{\left(\frac{B_z}{0.1\text{G}}\right)}^{-4/9}{\left(\frac{M_{\text{d}}+M_{\star }}{0.1M_{\odot }}\right)}^{1/3},$$(10)

where δ is the ratio of the initial density profile and the singular isothermal sphere (SIS, Shu 1977), and M_d is the disk mass. By comparing our density profile to the SIS, we take δ = 10, in agreement with Hennebelle et al. (2011), and the mean magnetic field strength within the disk as a proxy for the component B_z. For η_AD we take the azimuthally averaged value at the cell located farther away from the sink, within the disk selection.

The disk sizes agree within a factor of ≈2–3 with the prediction above (right panel of Fig. 7). We find roughly similar disk radii in all runs (Fig. 7), in agreement with this model predicting that the disk radius does not explicitly depend on the amount of angular momentum available (incontrast to the hydro case). The disk around primary sinks therefore appears to be set by magnetic regulation.

The deviation from the analytical prediction is similar to what has been found in C21: The disk radius is slightly overestimated analytically for runs μ = 2 and μ = 5, and it is underestimated when rotational support is significant (our SUPA and SUPAS runs, and run MU5ADf in C21 with μ = 5 and E_rot ∕E_grav = 5%). This underestimation may be due to the spiral patterns in disks when support from rotation or turbulence is strong, as the disk becomes gravitationally unstable.

Fig. 8 Column density maps at t ≈ 50 kyr. The integration was made in the direction of the angular momentum vector. From left to right: runs NOTURB, SUPA, SUPAS, and SUBA. White dots represent sink particles.

Fig. 9 Azimuthal median of the column density as a function of the radius, for the four runs, at t = 50 kyr. Coloured circles indicate the radius derived from our disk definition.

3.3.3 Column density maps

Our disk definition is physically motivated (see Sect. 2.4) and allows a comparison with previous numerical studies. Nevertheless, even though Keplerian disks can be traced from the position-velocity diagram (e.g. Cesaroni et al. 2005), the other criteria can hardly be verified observationally. For this purpose, Fig. 8 shows the column density maps at t = 50 kyr. The column density of individual disks is higher than that of their surroundings by at least one order of magnitude. Spiral arms are visible, as is circumbinary gas. Figure 9 shows the azimuthal median of the column density as a function of the radius. Colored circles indicate the radius that was determined using our disk definition. We chose a median rather than a mean in order to reduce the effect of non-axisymetries due to fragments or sinks. A smooth but visible transition (in slope and values) is observed around the disk radius we derived. For runs SUPA and SUPAS, the dense circumbinary gas makes determining the primary disk more difficult.

3.3.4 Disk fragmentation

As shown above, the origin of sink formation is disk fragmentation. This occurs in high-density gas at the extremity of spiral arms, similarly to the nonmagnetic case studied in Mignon-Risse et al. (2020). Spiral arms are found in all runs and at (nearly) all times. We determined the Toomre parameter (see Kratter & Lodato 2016 for a review of disk fragmentation) value before the development of spiral arms because it indicates the unstable state to axisymmetric perturbations that lead to spiral arm development (as is often done in the literature, see, e.g., Kratter & Matzner 2006; Klassen et al. 2016; Ahmadi et al. 2018). We computed it as in Mignon-Risse et al. (2020) and do not display it here for conciseness. Disks are Toomre-unstable with typically Q ≲ 0.7, except in their outer parts (r ≳ 100 AU). As shown inFigs. 2 and 12, disks are (largely) dominated by thermal pressure rather than magnetic pressure (β > 1). This means that it is sufficient to compute the thermal Toomre (without the magnetic pressure).

However, spiral arms do not indicate fragmentation or the formation of a multiple stellar system. In run NOTURB, no fragment forms. In run SUBA, the first fragments form at a sink age of about ~ 32 kyr, but fall and merge back onto the primary disk. In run SUPA, the first fragment formation occurs at a sink age ~ 18 kyr, and in run SUPAS, it occurs before 16 kyr. These differences suggest that the interplay of turbulence and magnetic fields may affect disk fragmentation and sink formation.

Several criteria are frequently used in the literature to address the origin of disk fragmentation in addition to the Toomre Q parameter. Klassen et al. (2016) used the Gammie criterion (Gammie 2001), aiming at comparing the cooling time to the orbital time. Even though we find that the local cooling time is generally shorter than the orbital time (even in our nonfragmenting run NOTURB, in agreement with Klassen et al. 2016 and using the same procedure), radial radiative flux is propagating in the disk to heat it up at a similar rate. This means that cooling is not responsible for fragmentation here. The Gammie model is well adapted to disks possessing local cooling processes, while our disks heat and cool radiatively. Moreover, Gammie (2001) assumed that the disk was cool and very thin, while we obtain an aspect ratio of typically 0.1–0.2. Furthermore, his local model is only applicable if c_s ∕(rΩ) ≪ 0.12 (Eq. (26) of Gammie 2001), where Ω is the orbital frequency. With c_s ≈ 10⁵ cm s⁻¹, r ≈ 100 AU and Ω ~ 10⁻¹¹ s⁻¹, we obtain c_s∕(rΩ) ~ 10 ≫ 0.12. We therefore argue that the model of Gammie (2001) does not apply to the massive protostellar disks formed in this study.

We investigated spiral arm collision as a possibility for forming fragments. Figure 10 shows the gas density as a function of position (in a fixed direction in the disk plane) and time in run SUPA. Spiral arms are visible at all times as diagonal lines of enhanced density in the [x, t] plane. A spiral arm collision event is indicated by the two horizontal lines (occurring in a similar fashion as in Oliva & Kuiper 2020). It creates a region of enhanced density, that is, a fragment, and decouples it from its parent arm. Collisions are observed in runs SUBA and SUPAS as well. This process is favored when the spiral arms can extend sufficiently far away from the primary star radially, which does not occur in run NOTURB. The rapid growth of the spiral arms in the turbulent runs is shown in the left panel of Fig. 7. When turbulence is included, it adds angular momentum for the spiral arms to extend (Hennebelle et al. 2016b, 2017). The initial distribution of angular momentum computed with respect to the center of the domain increases with the distance. This means that the less turbulent the core, the longer it takes for gas with a similar angular momentum to reach the center. This explains why fragments form earlier in run SUPAS than in runs SUPA and SUBA (see below a comparison between runs SUPA and SUBA). To conclude, we find that disk fragmentation is due to spiral arm collision and favored by turbulence.

In order to understand the collapse of a fragment, we followed the properties of the fragment in SUPA presented above. Cells with a density higher than 10⁻¹²g cm⁻³ and a distance to the primary sink larger than 250 AU (to avoid the disk) were selected as part of the fragment. Figure 11 shows the mean density, mean sound speed, size, and ratio of the free-fall time to the thermal sound-crossing time of the fragment from its formation time (t =46.05 kyr) to the sink formation time (t = 46.71 kyr). The free-falltime was computed as in Eq. (4) with the fragment mean density, and the sound-crossing time was equal to 2R∕c_s, with 2R the estimate for the fragment diameter and c_s its mean sound speed. This ratio gives an order-of-magnitude estimate of the fragment stability to density perturbations: a ratio lower than one indicates that it is gravitationally unstable. This estimate only accounts for thermal pressure because it largely dominates magnetic and radiative pressures locally. While the sound speed and density are found to be roughly constant with time, the fragment radius increases by a factor ~ 2.5 compared to its initial value. The fragment has accreted mass. Moreover, these quantities directly dictate the evolution of the ratio of the free-fall time and sound-crossing time. It is proportional to $c_{\text{s}}/(\sqrt{\rho }R)$, thus it tends toward an unstable state as R increases. Figure 11 shows that it is slightly unstable during two distinct epochs. The sink forms during the second epoch, indicating that part of the fragment has become bound, in addition to being Jeans unstable.

The previous analysis of fragment formation does not allow us to distinguish runs SUBA and SUPA, that is, whether there is an effect from magnetic fields on disk fragmentation. The fragments formed in SUBA merge back onto the disk and SUPA produces a companion, while both runs have the same level of turbulence and therefore of angular momentum. Moreover, the first fragment forms significantly later in run SUBA than in SUPA. These two observations mean that the magnetic braking is stronger in SUBA than in SUPA (see Fig. 4 and Sect. 3.6). It takes longer for larger-scale gas – which carries moreangular momentum – to fall onto the disk and spiral arms, delaying the first fragment formation. Similarly, magnetic braking slows down the gas rotation and makes the fragment formed in SUBA fall onto the disk.

To sum up, turbulence brings additional angular momentum to the ubiquitous spiral arms, which do not grow significantly in the non-turbulent run. This additional angular momentum favors their radial extension, their subsequent collision with another spiral arm and therefore the formation of fragments. When turbulence is sub-Alfvénic, magnetic braking delays this process and drives the inward migration of fragments, preventing long-lived stellar systems to emerge.

Fig. 10 Modified logarithmic color map of the density as a function of the position and time in a given direction in the disk plane in run SUPA. Spiral arms appear as diagonal lines of enhanced density. We are particularly interested in the spiral arm collision at t = 46.05 kyr, indicated by the two horizontal lines located at ±0.04 kyr from the collision time.

Fig. 11 Properties of the fragment leading to secondary sink formation (t = 46.71 kyr) in run SUPA. The density, radius, and sound speed are normalized to the value at the formation time. All quantities are averaged over the fragment volume.

3.3.5 Characteristic velocities and magnetic field components

Figure 12 shows the radial profile of the azimuthally averaged characteristic velocities in the disk selection, using the disk plane as the (u_r, u_θ) plane. Overall, the azimuthal velocity agrees with a Keplerian profile. It is slightly super-Keplerian in turbulent runs (SUPA, SUPAS, and SUBA) at radii ≳60 AU. In all runs, it becomes sub-Keplerian as the radius decreases. The reason is that the gravitational field is dominated by the central object and diminished by the sink softening length mentioned in Sect. 2.3. The disks are roughly Keplerian, therefore the infall velocity is far lower than the free-fall velocity and typically lower than 1 km s⁻¹. The rotation motions and infall motions beyond the disk are supersonic. The cells in the primary disk plane in the binary systems appear to be close to Keplerianity up to ~ 1000 AU (not shown here) when measured at t = 50 kyr (i.e., when the secondary sink is much less massive than the primary sink). In the absence of strong stellar activity from one component, they might be identified as large disks (Johnston et al. 2015).

As shown in the last row of Fig. 2, all primary disks have plasma beta β >1 (thermally dominated). In density maps and Fig. 12 we observe that the disk radius is close to the point atwhich a change from a thermally dominated to magnetically dominated region is observed. We argue that a physically motivated criterion for the identification of individual disks is the plasma beta with β ≳1 (equivalent to P_th ≳ P_mag) because $c_{\text{s}}/v_{\text{A}}=\sqrt{\gamma \beta /2}~\sqrt{\beta }$, especially atlate times. In our simulations, this means that the disk size is regulated by ambipolar diffusion, in contrast to the ideal MHD case (Masson et al. 2016, C21). This criterion alone is not sufficient, however, because of the thermally dominated (β > 1) filaments and parts of the circumbinary disk (see Sect. 3.5), but it works well in the vicinity of the stellar object. The filaments can be discarded by an additional criterion based on rotation.

Figure 13 displays the azimuthally averaged magnetic field components using the same coordinates as in Fig. 12. We selected cells in the disk plane, but these were not restricted to the disk selection in order to probe the outer regions as well. Strikingly, the vertical component B_z dominates for r ≲ 50 AU in runs NOTURB and SUBA, and in run SUPA most of the computational time. A dominantly poloidal magnetic field is a necessary condition to launch centrifugal jets (Blandford & Payne 1982). In runs SUPA and SUPAS, the magnetic field components show more variations with respect to time than in runs NOTURB and SUBA. We observe many occurrences in which B_z and B_r have opposite evolutions. Changes between B_z and B_r could be explained by the orbital motion of the primary sink, which occurs at super-Alfvénic velocities. A change from B_z to B_r can be attributed to magnetic field lines lagging behind the sink, while ambipolar diffusion would favor a return to B_z, as in runs NOTURB and SUBA. Following the method from C21, we defined the orbital Elsasser number for ambipolar diffusion Am as the ratio of the orbital time τ_orb and the ion-neutral collision time $t_{\text{in}}={\eta }_{\text{AD}}{\text{c}}^2/(4\pi v_{\text{A}}^2),$ $$\text{Am}=\frac{4\pi v_{\text{A}}^2{\tau }_{\text{orb}}}{{\text{c}}^2{\eta }_{\text{AD}}}.$$(11)

Taking the values in the vicinity of the primary sink, we have v_A ~ 0.1 km s⁻¹ (see Fig. 12), η_AD = 0.1, and τ_orb ≃ 1 kyr (Fig. 7), which gives Am ~ 0.4. This is of the order of unity (we recall that it is highly dependent on v_A, which increases away from the sink) and indicates that, indeed, kinematical effects compete with ambipolar diffusion.

At larger radii, including within the disk radius, the toroidal component B_ϕ dominates in all runs because the magnetic field lines are twisted by the disk rotation. Eventually, the radial component B_r dominates at even larger radii. It has been produced by the magnetized, collapsing gas (and the streamers), pulling the field lines that fan out at infinity and form an hourglass shape (e.g., Galli & Shu 1993).

Fig. 12 Azimuthally-averaged radial and azimuthal velocities, Alfvén speed, isothermal sound speed, free-fall velocity and Keplerian velocity as a function of the radius at t = 50 kyr in the disk selection. The vertical line indicates the disk radius plotted in Fig. 7 (which encapsulates 90% of the diskmass). Top row: run NOTURB (left), SUPA (right). Bottom row: run SUPAS (left), SUBA (right). Discontinuities are due to the disk selection being on a cell-by-cell basis, without connectivity criterion.

3.4 Secondary disk (runs SUPA and SUPAS)

We focused on the disk around secondary sinks in run SUPA and SUPAS in comparison with the primary disk in the same run (Sect. 3.3). First, they rotate in the same direction as the primary disks. Figure 14 shows the radial profile of the azimuthally averaged characteristic velocities within the cells corresponding to our disk selection criteria applied to the secondary sink environment. Similarly to the primary disks, their rotation profile is consistent with Keplerian rotation and they have c_s > v_A. Their plasma beta is shown in Fig. 15 (which is a slice centered on the center of mass of the two sinks) and shows how similar the two disks are. The region in which the azimuthal velocity v_ϕ is inverted corresponds to the closest part of the primary disk that dominates the azimuthal average (seen also around the primary sinks at later times than displayed in Fig. 12). This also gives an upper limit to the secondary disk radius, which complements the transition radius at which β ~ 1.

In both runs, the toroidal component of the magnetic field dominates most of the secondary disk region. At small radii (< 50 AU), the dominant component is not constant with time, as seen around the primary sinks (Sect. 3.3.5). The coevolution of B_z and B_r likely results from the same mechanism, that is, a B_z component increased by ambipolar diffusion, but turning into B_r as the secondary sink moves at a super-Alfvénic speed. Nonetheless, the temporal evolution does not show a clear pattern that would permit us to link the observations described above with characteristic timescales (e.g., the orbital period). By the end of the run SUPA, the secondary disk becomes dominated by B_z most of the time, similarly to the primary disks in runs NOTURB, SUPA, and SUBA. While the long-term evolution of this system is beyond the scope of this study, we may expect a magnetic structure that is favorable to outflow launching to build more rapidly around the secondary sink in run SUPA than in run SUPAS.

The density maps (not shown here) and velocity profiles show (e.g., taking the radius at which v_ϕ drops as an upper limit, in Fig. 14), the secondary disk radius is found to be ~ 100–150 AU, in agreement with the transition radius β ~ 1. This is the same order of magnitude as was found for the primary disks, and it is consistent within a factor of 2 with magnetic regulation (Eq. (10) predicts ~200 AU in these cases).

Fig. 13 Azimuthally-averaged magnetic field components as a function of the radius at t = 50 kyr. Top row: run NOTURB (left), SUPA (right). Bottom row: run SUPAS (left), SUBA (right).

Fig. 14 Same as Fig. 14 for the secondary sinks in runs SUPA (left) and SUPAS (right) when their mass is 5 M⊙.

3.5 Circumbinary disk (runs SUPA and SUPAS)

Runs SUPA and SUPAS show the formation of binary systems. The primary disks are embedded in a rapidly evolving circumbinary disk-like structure.

We investigated whether the binary separation was consistent with hydrodynamical disk radii, whose size was set by the centrifugal barrier from initial angular momentum conservation, rather than magnetic regulation. Nonetheless, this calculation led to a constant value, while we observe elliptic orbits (with a factor of ~ 2 between the periastron radius and the apastron radius) whose distance increases with time, once integrated over several orbits. The core rotational energy is $E_{\text{rot}}=M_{\text{c}}{(R_{\text{c}}\Omega )}^2/2=M_{\text{c}}{(J/M_{\text{c}})}^2/(2R_{\text{c}}^2)$, where R_c is the radius and $J=\Vert {\displaystyle {\int }_{r<R_{\text{c}}}r}\times \rho v\text{d}V\Vert$ its angular momentum, centered on the primary sink. Equalling the rotational energy and the gravitational energy 3G M_c∕(5R) (assuming a uniform density), we obtain the hydro disk radius $$r_{\text{d},\text{hy}}\simeq 100\text{AU}\times {\left(\frac{J/M_{\text{c}}}{5\times {10}^{21}{\text{cm}}^{\text{2}}{\text{s}}^{\text{-1}}}\right)}^2{\left(\frac{M_{\text{c}}}{100M_{\odot }}\right)}^{-1}.$$(12)

The quantity J∕M_c is displayed in Fig. 16 as a function of the integration radius. We took the value for R = R_c. For run SUPA, J∕M_c ≃ 8 × 10²¹ cm² s⁻¹, hence r_d,hy ≃ 300 AU, and J∕M_c ≃ 10²² cm² s⁻¹ so r_d,hy ≃ 400 AU for run SUPAS. These values roughly meet the binary separation (see Sect. 3.3.1 and Fig. 7). Hence, the binary separation appears to depend on the initial turbulent velocity field. To gain generality, this experience should be repeated with other realizations of the initial turbulence, but this is beyond the scope of this study.

The circumbinary disk that surrounds the two sink + disk systems is about twice larger than the binary separation, and it evolves with time as the two disks interact and the second disk grows. At t = 50 kyr, which is close to the birth epoch of the two secondary sinks, the circumbinary disk mostly has β ≥ 1 in run SUPA (Fig. 2). As shown in Fig. 15, when the secondary sink mass is 5 M_⊙, most of the surrounding gas has become magnetically supported in run SUPA, while thermally supported gas infall continues actively in run SUPAS. The circumbinary disk does not appear to be an isolated structure from the two individual disks, but the fate of this accreting system (individual disks and circumbinary disk) deserves dedicated studies. In both runs, the binary system is surrounded by a mostly toroidal magnetic field, similarly to unitary systems (runs NOTURB and SUPA), as displayed in Fig. 13. Hence, the magnetic field geometry cannot be used here to distinguish between a unitary and a binary system.

Fig. 15 Plasma beta slices of 2500 AU in the rotation plane, centered on the center of mass, when the secondary sink mass is 5 M⊙. The magnetic field topology is overplotted. White corresponds to a dominantly toroidal field. Left panel: run SUPA. Right panel: run SUPAS.

3.6 Primary disk orientation

One objective of this work is to make progress on the question of whether disks and outflows align with core-scale magnetic fields. The disk normal might be expected to be preferentially perpendicular to the magnetic field because then the projection of the magnetic braking is smaller and it cannot fully prevent disk formation (Joos et al. 2012). However, this picture may change when turbulence and nonideal MHD are accounted for. Because both effects individually (see, e.g., Joos et al. 2013; Hennebelle et al. 2016a) and together (Lam et al. 2019) solve the magnetic catastrophe, we would expect the disk-magnetic field orientation to tend toward a random distribution. We investigated the specific angular momentum components, and the alignment between this vector and the large-scale magnetic field (along the x-axis) as a function of the Mach number and the magnetic field strength. The angular momentum vector computed on disk scales (< 10³ AU) reveals the disk orientation.

Figure 16 shows the specific angular momentum j defined as $j=J/M=\frac{1}{M}{\displaystyle {\int }_{r<R}r}\times \rho v\text{d}V$ as a function of the spatial scale R for three epochs: t = 0, 30, and 50 kyr. We took r as the radial vector with respect to the primary sink position, except for the first snapshot, where there is no sink and we took the center of the box. We recall that each run initially had a tiny rotational support (1%) of solid-body rotation aligned withthe x-axis. In our reference case, run NOTURB, the specific angular momentum is initially aligned with the magnetic field axis and remains so (within less than 6°, not shown here for readability). Figure 16 shows the increase in local specific angular momentum in the central part of the domain (while the total angular momentum remains constant) as collapse occurs due to angular momentum transport during the infall. The angular momentum set by the initial turbulence is dominated by its y-component. The dominating component of j (accounting for both rotation and turbulence) varies with the sphere radius over which it is computed as a consequence of the initial turbulent velocity field. The initial rotation, aligned with the x-axis, dominates at large scales (>10⁴ AU) in runs SUPA and SUBA (left and right panels), but not in run SUPAS, where it is lower. This means that the turbulent gas is in counter-rotation with respect to the initial solid-body rotation imposed. In the two runs with super-Alfvénic turbulence, SUPA and SUPAS (left and central panels of Fig. 16), the dominating components at disk scales (< 10³ AU) at t = 50 kyr are the initial dominating components at slightly larger scales. These are the x- and y-components in the subsonic run SUPA and the z- and y-components in the supersonic run SUPAS. Hence, the disk orientation is set by the initial angular momentum rather than by magnetic fields.

We focused on the effect of magnetic fields. The left and right panels of Fig. 16 only differ by the magnetic field strength: μ = 5 in run SUPA (left) and μ = 2 in run SUBA (right). At small scales, the component j_x in run SUPA is ≈2 times larger than in run SUBA. This is a consequence of the magnetic braking, that is, the transport of angular momentum outward. It prevents disk formation perpendicular to the magnetic field and favors configurations in which the angular momentum is misaligned with the magnetic field.

In Fig. 17, we show the angle between j (the disk normal) and the x-axis, which corresponds to the direction of the large-scale magnetic field. The orientation of the angular momentum varies significantly with the scale considered and with time. Our results for sub- and supersonic turbulences are similar to those of Joos et al. (2013), namely a strong misalignment between j and B. The orientation on small scales converges in time as the disk forms and increases in size. On larger scales, the orientation does not vary, except in the most turbulent run, SUPAS (Fig. 17). This change is likely due to the change in velocity field configuration, as it becomes dominated by the gravitational pull exerted by the center of mass. Comparing left and right plots of Fig. 17 shows that increasing the magnetic field strength up to the point at which the initial turbulence is sub-Alfvénic does not favor the alignment between j and B.

Overall, the disk normal in our simulations is misaligned (50–85 °, Fig. 17) with the large-scale magnetic field, largely because of the initial turbulence. If the disk formation were a large-scale process, we would expect the disk normal to align with the core-scale angular momentum. However, as shown in the middle and right panels of Fig. 16, j_y < j_z (to be distinguished from j_x, which is more affected by magnetic braking) at the disk scale, while j_y > j_z at core scales. The disk orientation here does not appear to be set by the angular momentum direction at core scales because the gas on core scales has not yet reached the center of the domain within our simulated time (≲ 80 kyr), whereas disk formation occurs within a few 10 kyr. This would indicate that disk formation is a “local” process, in agreement with the recent observations in the low-mass regime (Gaudel et al. 2020), but this needs first to be confirmed for other initial density profiles. Moreover, disk evolution (and orientation) may be affected by this core-scale angular momentum, but longer-time integration is required.

Fig. 16 Specific angular momentum $j=J/M=\frac{1}{M}$ ∫r<Rr × ρv dV as a functionof the sphere radius R for runs SUPA (left), SUPAS (middle), and SUBA (right). Components of the specific angular momentum vector are displayed at t = 0 kyr (dotted lines), t = 30 kyr (dashed lines), and t = 50 kyr (full lines). Time t = 0 kyr describes the initial conditions, t = 30 kyr is roughly the first sink formation epoch (a rotating structure is already present), and t = 50 kyr corresponds to a massive protostar surrounded by its accretion disk. Components along x, y, and z are shown in blue, orange, and green, respectively.

Fig. 17 Angle of the specific angular momentum j and the x-axis (the initial magnetic field orientation) as a function of the sphere radius R for runs SUPA (left), SUPAS (middle), and SUBA (right). The same time steps as Fig. 16 are pictured here. The angle is displayed at t = 0 kyr (dotted lines), t = 30 kyr (dashed lines), and t = 50 kyr (full lines).

4 Discussion

4.1 Comparison with previous works

As mentioned previously, this work extends the study of C21 to a turbulent medium and with a hybrid radiative transfer method. In our nonturbulent run NOTURB, we have obtained a final sink mass of M ≃ 15.8 M_⊙. This value is very similar to what has been found with the FLD method in C21 (their nonideal MHD run with μ = 5 and E_rot ∕E_grav = 1%). The disk radius obtained in run NOTURB also compares well with C21, with ~ 100 AU. It shows that in a magnetized environment, the radiative feedback method does not seem to regulate the sink mass or the disk radius, up to M ≃ 15.8 M_⊙. The effect of magnetic fields and ambipolar diffusion is undeniable, however. With high-resolution simulations including ohmic dissipation (and no ambipolar diffusion), Kölligan & Kuiper (2018) obtained disks of up to 10³ AU. While the Hall effect is expected to reduce or increase the disk size depending on the alignment of the rotation axis and the magnetic field axis, this suggests that the strong regulating role of ambipolar diffusion (Hennebelle et al. 2016a) is unique among nonideal MHD effects. Taking advantage of this role, we propose that a plasma beta β > 1 (thermal pressure exceeding magnetic pressure) is a good indicator for distinguishing individual disks in the vicinity of protostars.

We report the presence of accretion streamers (already observed in the low-mass regime, see Pineda et al. 2020), which are reminiscent of the accretion channels found by Seifried et al. (2015) and the bridges in the study of Kuffmeier et al. (2019) and Kuffmeier et al. (2020). In agreement with Seifried et al. (2015), ram pressure dominates magnetic pressure, but we further show that magnetic pressure is also dominated by thermal pressure. These structures develop even in the nonturbulent run NOTURB, and they seem to be associated with magnetic field effects and are not a pure consequence from turbulence, as shown in Kuffmeier et al. (2019). We note that the width of the accretion streamers appears to depend on the initial turbulent level, however: the stronger the turbulence, the wider the streamer. This effect could be a consequence of magnetic turbulent reconnection that occurs in the envelope and increases with the turbulent level (e.g., Santos-Lima et al. 2013). The physical origin of these accretion streamers should be investigated in dedicated studies.

The present work confirms that disk formation does not occur preferentially perpendicular to the core-scale magnetic fields, but its orientation is likely driven by turbulence (even sub-Alfvénic). This is in agreement with the study of Machida et al. (2019) in the low-mass regime, in which they varied the angle between the rotation axis and the magnetic field direction. They observed that the disk plane is mainly set by their initial rotation (i.e., specific angular momentum) axis, even with an initially strong magnetic field (μ = 1.2).

Very few works on massive star formation have reported the formation of a binary system. Here, we obtained mass ratios of ≈ 1.1–1.6 between the two sinks. Balanced mass ratios like this (of the order of unity) have been obtained in the radiation-hydrodynamical simulations of Krumholz et al. (2009), where they could be integrated for longer times, with final masses of 41.5 M_⊙ and 29.2 M_⊙ and a separation of 1590 AU. For comparison, the binary separation is 350–600 AU in run SUPA and 400–700 AU in run SUPAS (the orbits are elliptic, and the separation increases slightly with time). The studies of Meyer et al. (2018), focused on the formation of spectroscopic binaries, and Oliva & Kuiper (2020) on disk fragmentation are hardly comparable with the present work because they only used a sink particle algorithm for the primary star.

4.2 Comparison with observations

We first compare our findings in terms of disk radius with observations. The disk radius of HH80-81, estimated to be ~ 291 AU (Girart et al. 2018) or Orion source 1 with ≲50 AU (Matthews et al. 2010; Ginsburg et al. 2018), agree better with the magnetically regulated indivual disk radii we obtain than with purely hydrodynamical disks (see, e.g., Kuiper et al. 2011; Mignon-Risse et al. 2020). This is also in line with the upper limit of 125 AU set on the disk of S255IR SMA1 (Liu et al. 2020). The binary separation is linked to the centrifugal radius that can be derived from the initial turbulent field. IRAS 16547-4247 is a rare and recent case of a massive protostar binary. Recent measurements with ALMA indicate jets, a circumbinary disk of radius ~ 2500 AU, individual disks on a scale of ~100 AU, and a binary separation of 300 AU (Tanaka et al. 2020). These order-of-magnitude estimates are consistent with runs SUPA and SUBA, except that they observe indications of counter-rotating disks. This indicates preferentially another origin than disk fragmentation.

Furthermore, Aizawa et al. (2020) have studied the disk orientation in five nearby star-forming regions and observed a random orientation, consistent with the disk plane being set by turbulence. This agrees with our work. This aspect is also a preliminary step to identify the orientation of outflows (naively expected to be perpendicular to the disk), which has been observationally investigated (e.g., Hull et al. 2020).

4.3 Limitations

We now discuss some of the limitations of our approach. As many other studies of massive star formation (e.g., Krumholz et al. 2009; Kuiper et al. 2010a; Commerçon et al. 2011a), we chose a high-mass prestellar core for our initial conditions, consistent with the turbulent core model of McKee & Tan (2003). This approach allowed us to reach a finest resolution of 5 AU to capture the disk and the physical scales of interest (the Jeans length and the disk scale height) except very close to the sink, which is rarely done in large-scale simulations without zoom-in procedures. This condition is strengthened by the nonideal MHD frame, which leads to smaller disks than in the hydrodynamical case (see, e.g., Hennebelle et al. 2016a; Masson et al. 2016). As claimed by various models, such as the global hierarchical collapse model (Vázquez-Semadeni et al. 2009) and the inertial-inflow model (Padoan et al. 2020), and supported by various observations (Schneider et al. 2010), the large-scale dynamics likely plays a major role in the formation of massive stars, and the isolation we have assumed may be an oversimplification, in particular for the turbulence as it is modeled in this paper. We leave this to further work.

Protostellar heating could suppress disk fragmentation or promote it in the shielded disk midplane because heating increases the Jeans length and stabilizes the gas against gravitational collapse. However, most of the protostellar radiation is absorbed at the disk inner edge and does not directly heat the gas located farther away. Very high resolution is required to resolve the photon mean free path, derive the exact temperature structure, and conclude on the effect of protostellar heating on disk fragmentation. A cell of gas density 10⁻¹¹g cm⁻³ and opacity to ultraviolet radiation ~10²cm² g⁻¹ would require a spatial resolution ofless than 10⁻³ AU. However, this highly depends on the opacity, hence on the source frequency and on the dust density.

Our sink accretion criterion relies on the local Jeans density (as in C21). Investigating the effect of the accretion method (e.g., a density threshold) would be an asset, but this is beyond the scope of this paper.

We enforced gas-radiation decoupling within the sink volume for outflow physics purposes (see Paper II). The star exerts a stronger radiative force onto the gas in the first absorption region, that is, at the sink accretion radius (~ 20 AU). This possibly shifts the disk inner edge (i.e., the edge at 20 AU receives a direct stellar radiative force, perturbing the gravitation–centrifugal equilibrium; see Appendix A).

Finally, we used a constant dust-to-gas ratio (1%) within the simulated volume throughout this paper. Nonetheless, dust is the main contributor to the medium opacity. Grain growth and sedimentation are expected to occur, affecting this dust-to-gas ratio and therefore the opacities that couple the dust-gas mixture and radiation. Furthermore, dust grains are charge carriers, which means that the ionization degree would vary and the nonideal MHD resistivities with it. Dust dynamics should therefore be integrated in collapse calculations (Lebreuilly et al. 2019; Lebreuilly et al. 2020). This would allow one to obtain a dust size distribution that varies dynamically, affecting the dust-and-gas mixture temperature and the effect of magnetic fields, which in turn sets the radial equilibrium of the disk.

5 Conclusions

We have conducted four numerical simulations of massive prestellar core collapse including ambipolar diffusion and a hybrid radiative transfer method. This led to the formation of thermal pressure-dominated disks, rather than magnetic pressure dominated disks. We included an initial velocity field consistent with turbulence, and varied the Mach and Alfvénic Mach numbers to consider four runs without turbulence (NOTURB), super-Alfvénic-subsonic turbulence (SUPA), super-Alfvénic-supersonic turbulence (SUPAS), and sub-Alfvénic-subsonic turbulence (SUBA). We summarize our results below:

• 1.

  Even in the absence of turbulence, asymmetries are naturally caused by streamers (thermally dominated filaments slightly denser than their surroundings) that connect onto the disk off the disk plane, and by the interchange instability, which redistributes magnetic flux at the disk edge (see Appendix B), breaking the axisymmetry. Regarding the streamers, the seed is likely numerical in run NOTURB, but not in the non-axisymmetric runs SUPA, SUPAS, and SUBA, where they are present as well;

• 2.

  Keplerian disks form in all runs. They have typical radii of 100–200 AU around individual stars and are consistent with magnetic regulation. In the super-Alfvénic runs, they are located within a larger rotating structure (circumbinary disk, see below). In this case, the rotation profile is close to Keplerian rotation within a few hundred AU;

• 3.

  We report the formation of stable binary systems when turbulence is super-Alfvénic. They form from disk fragmentation at the extremity of spiral arms rather than initial (core) fragmentation, and follow spiral arm collision. Their binary separation is between 300 AU and 700 AU and may be linked to the initial angular momentum (i.e., amount of rotation) carried by the turbulent velocity field;

• 4.

  We have assessed the misalignment between the disks and core-scale magnetic fields. The disk orientation appears to be set by the initial angular momentum at scales ≲10⁴ AU only, in agreement with the previous numerical study of Machida et al. (2019). The streamers are located in a plane perpendicular to the magnetic field in all runs, but do not affect the disk formation process;

• 5.

  A plasma beta (ratio of the thermal pressure and magnetic pressure) higher than one indicates structures of interest such as the individual disks and infalling streamers.

We presented disk accretion as the only accretion mechanism for massive protostars, in contrast to alternative mechanisms such as filamentary accretion or stellar mergers (Bonnell et al. 2001). The case of radiative Rayleigh-Taylor instabilities was not addressed here because of the debate on the required resolution to obtain them, but we refer to Krumholz et al. (2009), Kuiper et al. (2012), Rosen et al. (2016), and Mignon-Risse et al. (2020).

Our work confirms that multiplicity may be linked to medium turbulence. Depending on the models, this turbulence is expected to be higher in massive star-forming regions (due to radiative outflows and photoionization from other stars, and inflow from large scales), which may produce a higher stellar multiplicity than for low-mass stars.

Acknowledgements

This work was supported by the CNRS “Programme National de Physique Stellaire” (PNPS). The numerical simulations presented here were run on the CEA machine Alfvén and using HPC resources from GENCI-CINES (Grant A0080407247). The visualisation of RAMSES data has been executed with the OSYRIS python package. RMR particularly thanks T. Foglizzo for his insights on the interchange instability.

Appendix A Luminosity injection into the sink particle volume: disk size

In this appendix, we investigate the effect on the disk size of the radiative transfer method and of the kernel function when the luminosity is deposited within the sink volume. This is mainly motivated by the need to properly model photon escape within the sink volume to study radiative outflows (see Paper II).

Fig. A.1 Evolution of the disk radius as a function of the sink age. We varied the radiative transfer and the luminosity injection methods.

The simulations are similar to run NOTURB presented in the main text, therefore they include nonideal MHD (ambipolar diffusion), but no turbulence. We ran four simulations: two with the flux-limited diffusion (FLD), and two with the hybrid radiative transfer approach (HY). For each method, we tested two injection kernels: either the luminosity was deposited uniformly over the sink volume (uniform), or only over the central oct (peaked).

As shown in Fig. A.1, we obtain similar disk sizes in all runs, except in run HY uniform, where the disk radius is larger by ~20 AU. This corresponds to the sink accretion (and luminosity injection) radius. When the luminosity is deposited uniformly up to a radius of 20 AU, this leads to an additional repulsive force (the direct stellar radiative force) that is exerted onto the gas over a radius equal to the sink luminosity injection radius, which is 20 AU. Hence, a uniform luminosity injection with the hybrid method likely shifts the disk toward a slightly larger radius.

Nevertheless, a uniform injection of luminosity within the sink volume is not physically satisfying. The M1 radiative flux that powers the radiative force indirectly depends on the local radiative energy gradient. If the injection is uniform over the sink volume, radiative energy is more absorbed in the central cells (which are located in dense gas) than above and below the disk plane (where lower-density gas is located). This results in a radiative flux oriented toward the central cells and consequently in a spurious radiative force oriented toward the central cells, from above and below the disk plane. For this reason, we did not adopt a uniform luminosity injection function in this paper, but rather set the sink volume as entirely optically thin.

Appendix B Interchange instability

Figure B.1 shows that a pocket of magnetized plasma is released from the disk edge. This occurs several times in the simulations, but is more difficult to distinguish in the turbulent runs. In this section we determine whether the interchange instability (also called magnetic Rayleigh–Taylor instability), which is a convective instability that redistributes the magnetic flux, is responsible for this.

The instability occurs in the y−direction if (see, e.g., Lovelace & Scott 1981; Kaisig et al. 1992) $$\frac{\Vert \frac{\partial ((\gamma -1)P+B_{\text{x}}^2/2)}{\partial y}\Vert }{\left((\gamma -1)c_{\text{s}}^2+v_{\text{A}}^2\right)\Vert \frac{\partial \rho }{\partial y}\Vert }>1,$$(B.1)

where x is the normal direction to the disk, and $v_{\text{A}}=B/\sqrt{\rho }$ is the Alfvén velocity. The condition of instability is roughly given by the balance between the density gradient set by gravity and the (total) pressure gradient. We derived the growth rate ω =Im(N) analogously to the Brunt-Väisäla frequency (which is a frequency associated with convective instabilities) from $$N^2=\frac{1}{1+\alpha }\left(\frac{\gamma -1}{\gamma }\frac{\partial s}{\partial y}+\alpha \frac{\partial \log (B/\rho )}{\partial y}\right)g_{\text{eff}},$$(B.2)

where we defined that $\alpha \equiv v_{\text{A}}^2/c_{\text{s}}^2$, $s=\frac{1}{\gamma -1}\ln \left(P{\rho }^{-\gamma }\right)$ is the normalized gas entropy, and $g_{\text{eff}}=g-v_{\varphi }^2/r$ is the effective gravity at radius r (sum of the gravitational and centrifugal accelerations).

The bottom panel of Fig. B.1 shows the square root of N² in the y direction in the disk plane, varying the x and z coordinates of the origin among the four closest cells to the sink center. Zones where this value is purely imaginary are unstable, which correspond to the disk edge in multiple directions. The growth rate at the disk edge is ω ≈ 70 kyr⁻¹ for the origin cells of coordinates [3, 2] AU (origin 1) and [−2, 2] AU (origin 2), so that the timescale for the instability to develop is τ_instab ≈ 14 yr. The unstable zone is ≈20 AU wide, in which the gas is flowing at a radial velocity v_r ≈ 0.8 km s⁻¹, so that the advection timescale is τ_adv ≈ 120 yr. Because τ_instab ≲ τ_adv∕3 (Foglizzo et al. 2006), this is consistent with the interchange instability being at work. When taking the cells [3, −3] (origin 3) and [−2, −3] AU (origin 4) as the origin of the profile, it is less clear whether this part of the disk edge is stable. Hence, the small unstable part of the disk edge may explain why this instability is less visible than in Krasnopolsky et al. (2012) and too faint to be observable, even thoughthe interchange instability is a good candidate.

Fig. B.1 Interchange instability in run NOTURB. Top left and right panels: plasma β at t = 45.25 kyr (before the instability develops) and t = 45.37 kyr. Bottom panel: square root ω of N2 (see Eq. (B.2)): Im(ω) gives the growth rate of the interchange instability. We compute it in the y-direction in the disk plane, at t = 45.25 kyr, taking the four closest points to the sink center as origin.

References

All Tables

All Figures

Fig. 1 Density slices perpendicular (first and second row, zoomed by a factor of 10) and parallel (third and fourth row, zoomed by a factor of 5) to the disk plane, at t ≈ 50 kyr. Streamlines corresponding to magnetic field lines and arrows corresponding to the velocity field are overplotted. Columns from left to right: Runs NOTURB, SUPA, SUPAS, and SUBA. The mass density ρ = 10−19 g cm−3 corresponds to the particle density n = 2.6 × 104 cm−3 and ρ = 10−11 g cm−3 to n = 2.6 × 1012 cm−3. White dots represent sink particles.

In the text

	Fig. 2 Same as Fig. 1, but with plasma beta (ratio of the thermal pressure and the magnetic pressure) slices. The gas moves along magnetic field lines until it forms thermally dominated infalling filaments. Disks have β > 1 as well.
In the text

	Fig. 3 Density-magnetic field strength histograms at t = 50 kyr. From left to right: run NOTURB and SUPA (top), SUPAS and SUBA (bottom).
In the text

	Fig. 4 Evolution of angular momentum transported by the outflows (left panel), magnetic fields (middle panel), and by gravitation (right panel) within a cylinder of radius and height 1000 AU.
In the text

	Fig. 5 Primary and secondary sink mass (left panel) and accretion rate (right panel) as a function of primary sink age for the four runs. The accretion rate is plotted for one sub-Alfvénic (NOTURB)and one super-Alfvénic (SUPA) run for readability. It has been smoothed over periods of ~ 1 kyr.
In the text

	Fig. 6 Primary sink mass against the disk mass (left panel) and disk mass as a function of time (right panel) for the four runs. In the right panel, colored circles indicate the secondary sink formation epoch.
In the text

	Fig. 7 Disk radius (left panel) and ratio of the disk radius and the theoretical value (right panel, Eq. (10)) as a function of time for the four runs. Colored circles indicate the secondary sink formation epoch.
In the text

	Fig. 8 Column density maps at t ≈ 50 kyr. The integration was made in the direction of the angular momentum vector. From left to right: runs NOTURB, SUPA, SUPAS, and SUBA. White dots represent sink particles.
In the text

	Fig. 9 Azimuthal median of the column density as a function of the radius, for the four runs, at t = 50 kyr. Coloured circles indicate the radius derived from our disk definition.
In the text

	Fig. 10 Modified logarithmic color map of the density as a function of the position and time in a given direction in the disk plane in run SUPA. Spiral arms appear as diagonal lines of enhanced density. We are particularly interested in the spiral arm collision at t = 46.05 kyr, indicated by the two horizontal lines located at ±0.04 kyr from the collision time.
In the text

	Fig. 11 Properties of the fragment leading to secondary sink formation (t = 46.71 kyr) in run SUPA. The density, radius, and sound speed are normalized to the value at the formation time. All quantities are averaged over the fragment volume.
In the text

Fig. 12 Azimuthally-averaged radial and azimuthal velocities, Alfvén speed, isothermal sound speed, free-fall velocity and Keplerian velocity as a function of the radius at t = 50 kyr in the disk selection. The vertical line indicates the disk radius plotted in Fig. 7 (which encapsulates 90% of the diskmass). Top row: run NOTURB (left), SUPA (right). Bottom row: run SUPAS (left), SUBA (right). Discontinuities are due to the disk selection being on a cell-by-cell basis, without connectivity criterion.

In the text

	Fig. 13 Azimuthally-averaged magnetic field components as a function of the radius at t = 50 kyr. Top row: run NOTURB (left), SUPA (right). Bottom row: run SUPAS (left), SUBA (right).
In the text

	Fig. 14 Same as Fig. 14 for the secondary sinks in runs SUPA (left) and SUPAS (right) when their mass is 5 M⊙.
In the text

	Fig. 15 Plasma beta slices of 2500 AU in the rotation plane, centered on the center of mass, when the secondary sink mass is 5 M⊙. The magnetic field topology is overplotted. White corresponds to a dominantly toroidal field. Left panel: run SUPA. Right panel: run SUPAS.
In the text

Fig. 16 Specific angular momentum $j=J/M=\frac{1}{M}$ ∫r<Rr × ρv dV as a functionof the sphere radius R for runs SUPA (left), SUPAS (middle), and SUBA (right). Components of the specific angular momentum vector are displayed at t = 0 kyr (dotted lines), t = 30 kyr (dashed lines), and t = 50 kyr (full lines). Time t = 0 kyr describes the initial conditions, t = 30 kyr is roughly the first sink formation epoch (a rotating structure is already present), and t = 50 kyr corresponds to a massive protostar surrounded by its accretion disk. Components along x, y, and z are shown in blue, orange, and green, respectively.

In the text

	Fig. 17 Angle of the specific angular momentum j and the x-axis (the initial magnetic field orientation) as a function of the sphere radius R for runs SUPA (left), SUPAS (middle), and SUBA (right). The same time steps as Fig. 16 are pictured here. The angle is displayed at t = 0 kyr (dotted lines), t = 30 kyr (dashed lines), and t = 50 kyr (full lines).
In the text

	Fig. A.1 Evolution of the disk radius as a function of the sink age. We varied the radiative transfer and the luminosity injection methods.
In the text

	Fig. B.1 Interchange instability in run NOTURB. Top left and right panels: plasma β at t = 45.25 kyr (before the instability develops) and t = 45.37 kyr. Bottom panel: square root ω of N2 (see Eq. (B.2)): Im(ω) gives the growth rate of the interchange instability. We compute it in the y-direction in the disk plane, at t = 45.25 kyr, taking the four closest points to the sink center as origin.
In the text