Collapse of turbulent massive cores with ambipolar diffusion and hybrid radiative transfer I. Accretion and multiplicity

(Abridged) Context. Massive stars form in magnetized and turbulent environments, and are often located in stellar clusters. Their accretion mechanism, as well as the origin of their system's stellar multiplicity are poorly understood. Aims. We study the influence of both magnetic fields and turbulence on the accretion mechanism of massive protostars and their multiplicity. Methods. We present a series of four Radiation-MHD simulations of the collapse of a massive magnetized, turbulent core of 100 $M_\odot$ with the AMR code Ramses, including a hybrid radiative transfer method for stellar irradiation and ambipolar diffusion. We vary the Mach and Alfvenic Mach numbers to probe sub- and superalfvenic turbulence as well as sub- and supersonic turbulence regimes. Results. Subalfvenic turbulence leads to single stellar systems while superalfvenic turbulence leads to binary formation from disk fragmentation following spiral arm collision, with mass ratios of 1.1-1.6. In those runs, infalling gas reaches the individual disks via a transient circumbinary structure. Magnetically-regulated, thermally-dominated (plasma beta $\beta>1$), Keplerian disks form in all runs, with sizes 100-200 AU and masses 1-8 $M_\odot$. The disks around primary and secondary sink particles share similar properties. We observe higher accretion rates onto the secondary stars than onto their primary star companion. The primary disk orientation is found to be set by the initial angular momentum carried by turbulence. Conclusions. Small (300 AU) massive protostellar disks as those frequently observed nowadays can only be reproduced so far in the presence of (moderate) magnetic fields with ambipolar diffusion, even in a turbulent medium. The interplay between magnetic fields and turbulence sets the multiplicity of stellar clusters. A plasma beta $\beta>1$ is a good indicator of streamers and disks.


Introduction
Massive stars (> 8 M ) are luminous, among the main sources of feedback at parsec and galactic scales, especially due to their explosion into supernovae (Larson 1974). Nonetheless, the conditions in which they form remain unclear. Indeed, this challenging problem offers observational issues: massive stars are rare, located far away from us ( 1 kpc) and in dense regions; and theoretical difficulties: in addition to magnetic fields and kinetic energy from turbulence (Tan et al. 2014), radiative energy is to be accounted for. At late stages, they are also expected to drive Hii region expansion (Keto 2007), bringing even more complexity to their environmental conditions of birth. Hence, while the accretion mechanism in low-mass star formation is becoming increasingly understood in terms of disk-mediated accretion, it is not settled yet for high-mass star formation. Furthermore, sev-eral outflow models rely on disk accretion, hence a first uncertainty on the origin of outflows follows from the uncertainty on 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, in particular whether it arises from core or disk fragmentation. These three questions: accretion, ejection, and fragmentation (as a cause for multiplicity), are tightly linked together and depend on common physical ingredients: magnetic fields, turbulence, and radiation, that need to be modelled 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 two last ingredients, namely turbulence and radiation, thus the first attempts to model massive star formation have built on this. First of all, in Article number, page 1 of 22 arXiv:2105.14543v1 [astro-ph.SR] 30 May 2021 A&A proofs: manuscript no. aanda the Competitive Accretion scenario (Bonnell et al. 2001), lowand high-mass protostars feed from a common gas reservoir after large-scale core fragmentation. In this case, the final stellar mass is uncorrelated 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 pre-stellar cores, stabilized against gravitational collapse by turbulent motions and magnetic fields. This model suffers from the rareness of high-mass pre-stellar cores , though candidates do exist (e.g., Nony et al. 2018), and may not be the most common procedure for massive star formation. 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. 2019) challenge this core-fed accretion scenario and give prominence to large-scale dynamics, either due to collapse (GHC) or to the inertial motions following supernova feedback (II). Those suggest the inclusion (and possible requirement) of turbulence. It is not sure yet whether accretion should be "clump-fed" or "core-fed" and occur via disks or turbulent filaments (see e.g., Rosen et al. 2019), but recent observations put (sparse but) increasingly convincing constraints on disk-mediated accretion, thanks to unprecedented angular resolution.
Current observational constraints on disks indicate radii of 20 to 330 AU (Patel et al. 2005, 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 can be explained by the protostar's environment (magnetic field strength, turbulence, geometry?). Disks can be subject to fragmentation, possibly leading in the formation of multiple stars gravitationally bound together. The massive hot core region G351.77-0.54 observed with ALMA at sub-40 AU resolution reveals twelve sub-structures within a few thousand AU with a broad range of core separations , consistent with thermal Jeans fragmentation of a dense core, and possibly with the Global Hierarchical Model (Vázquez-Semadeni et al. 2016). The aforementionned disk in HH 80-81 could be prone to fragmentation (Fernández-López et al. 2011) as well. There does not seem to be any general trend about the disk stability, but the advent of ALMA will increase the statistics. In the high-mass star-forming region IRAS 23033+5951, four mmsources are identified with the Northern Extended Millimeter Array (NOEMA) and the IRAM telescope, out of which two exhibit protostellar activity. Among those, one 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, hence clump accretion can be detected by luminosity outbursts. The process behind 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 seems, however, a recurrent mechanism in low-mass star formation, and relies on diskmediated accretion and star-disk interaction. Hence, it advocates the same accretion method for the formation of high-mass stars as well (Caratti o Garatti et al. 2017).
Altogether, despite the lack of systematic constraints, the presence of disks around young massive protostars (L < 10 5 L ) is now well-established (see the reviews by Beltrán & de Wit 2016, Beltrán 2020. Their 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, due to new polarimetric instruments. In a sample of 21 high-mass star-forming clumps, sub-parsec magnetic fields appear to be structured (Zhang et al. 2014). The hour-glass shape due to field lines being 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 µ = (M/φ)/(M/φ) crit is the mass-to-flux to critical mass-to-flux ratio, 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 , 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 1 (Chandrasekhar & Fermi 1953). Comparisons with magneto-hydrodynamical simulations have shown that fragmentation is consistent with turbulence dominating over 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 or to dominate over the turbulent energy (sub-alfvénic turbulence, Pillai et al. 2015) in several sources. Girart et al. (2013) have found equipartition between rotational energy, magnetic energy and turbulent energy in a fast rotating core, with µ = 6, indicating three mechanisms capable of slowing down the collapse.
Hence, observations agree on the presence of disk-like structures, with several occurrences of fragmentation and sizes of tens to hundreds AU. Magnetic fields have non-negligible strengths and their presence may affect both disk formation (and subsequently outflow launching), and its fragmentation into multiple stellar systems. Let us summarize the recent improvements made on 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 multi-dimensional 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 being radiated preferentially off the plane, ending in a small radiative force against the accretion flow, compared to the unidimensional view. Meanwhile, progress has been made in the low-mass star formation context with the inclusion of magnetic fields in numerical simulation in the ideal magneto-hydrodynamics (MHD) frame (e.g. Fromang et al. 2006). Many studies have shown that in a collapsing core, the flux-freezing condition leads to the accumulation of magnetic fields in the central region, inducing a strong magnetic braking and preventing disk formation (see e.g. Fromang 2008, andSeifried et al. 2011 in the high-mass regime). This is referred to as the magnetic catastrophe, since many disks are observed around low-and highmass protostars (Cesaroni et al. 2005). Three ingredients have been introduced and shown separately to allow for the formation of disks and reconcile numerical simulations and observations in that respect: misalignment between the rotation axis and the magnetic field axis, turbulence, and non-ideal MHD effects. Misalignment (Joos et al. 2012) and turbulence (Joos et al. 2013, Lam et al. 2019for low-mass stars, Seifried et al. 2012 for highmass) directly reduce the magnetic braking efficiency. Non-ideal (also called resistive) MHD effects, namely ambipolar diffusion (AD), Ohmic dissipation and the Hall effect provide a mechanism to limit the accumulation of magnetic fields strength and therefore the magnetic braking. AD is probably the most-studied non-ideal MHD effect, as it starts dominating at lower densities than the others, and indeed promotes disk formation in the low- ) and high-mass regime (Commerçon et al. 2021, hereafter C21). In several studies, non-ideal MHD appears as 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 on massive star formation have focused on the radiative transfer aspect, due to 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 the stellar radiation, also called irradiation, has been improved later on to track the higher-energy stellar photons and compute accurately the opacity during their first interaction 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. 2020. In the meantime, the common inclusion of radiative transfer and MHD in numerical codes has shown that both effects contribute to limit 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 modelling of magnetized disks requires non-ideal MHD effects to circumvent the so-called magnetic catastrophe (C21). In this work, we extend the C21 study and present the first numerical simulations including both a hybrid radiative transfer method and non-ideal MHD (namely, ambipolar diffusion), aiming at identifying the accretion conditions of massive protostars, as well as 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 non-idealized 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 on the primary, secondary and circumbinary disks properties and on the disk-magnetic field alignment. Our results and limitations are discussed in Sect. 4 and we conclude in Sect. 5.

Methods
In this section we present the set of equations that are 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 that has been done with the Flux-Limited Diffusion method for radiative transfer, but this time with an hybrid radiative transfer method and in a turbulent medium. An additional difference resides in the sublimation model we take here, and an optically-thin sink volume (see below).

Radiation magneto-hydrodynamical model
We integrate the equations of radiation-magneto-hydrodynamics (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 solve is 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 2 + 1/2B 2 + 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, P 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, P M1 is the M1 radiative pressure,Ė M1 is the stellar radiation injection term. Let us note that we inject radiative energy into the M1 module only with the primary sink, while other sinks (formed after the primary one) will radiate within the FLD module. The justification is twofold. First, in the M1 method the radiative fluxes from several sources sum up, whereas two radiation beams should not interact when crossing 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 one star's radiation is sufficient to do so.
The ambipolar electromotive force is equal to where η AD is the ambipolar diffusion resistivity. It depends on the density, temperature, and magnetic field strength. The resistivities are pre-computed using a chemical network to calculate the equilibrium abundances of the molecules (neutrals) and main charge carriers in conditions of pre-stellar core collapse (Marchand et al. 2016), depending on the density and temperature. Chemical equilibrium is assumed because the associated timescale is shorter than the free-fall time at the densities we consider.
The term κ P, ρcE M1 couples the M1 and the FLD methods via the equation of evolution of the internal energy We use the ideal gas relation for the internal specific energy ρ = C v T where C v is the heat capacity at constant volume.

Physical setup
We take similar initial conditions as C21. The free-fall time in the central plateau (which contains ∼15 M ) of the density profile is then where G is the gravitational constant andρ 0 is the density of the central plateau.
The initial core is threaded by a uniform magnetic field oriented along the x−axis. We set the magnetic field strength by the mass-to-flux to critical mass-to-flux ratio µ = (M/Φ) 0 /(M/Φ) crit where (M/Φ) 0 = M c /(πB 0 R 2 c ) and (M/Φ) crit = 0.53/(3π) √ 5/G (Mouschovias & Spitzer 1976). Strong (µ = 2, B 0 = 170 µG) and moderate (µ = 5, B 0 = 68 µG) magnetic fields are considered here. A drawback of this uniform distribution is that the mass-to-flux ratio decreases as µ ∼ 1/R and is larger in the inner parts of the core, with µ ≈ 50 in the central plateau (for runs with µ = 5) corresponding to a weakly magnetized medium. We expect, however, the central magnetic field strength to increase as B ∝ ρ 2/3 as the core contracts and before ambipolar diffusion starts dominating (at ρ∼10 −15 g cm −3 ). Thus, the magnetic field will play a dynamical role in the collapse (see C21).
An initial velocity dispersion is imposed to mimic a turbulent medium, and follows a Kolmogorov power spectrum P(k) ∝ k −5/3 , similar to Commerçon et al. (2011a). Phases are randomly sampled, and one realization is considered, i.e. we do not vary the velocity field but only its amplitude between two runs. The turbulence is not sustained but the turbulent crossing time (∼0.5 − 2 Myr with T = 20 K, see below) is significantly larger than the simulation time here, so it should not impact our results. A low level of (solid-body) rotation, E rot /E grav = 1%, is initially imposed around the x−axis, and dominates the specific angular momentum in subsonic runs. We consider four runs (see Table 1), varying the initial Mach number M and Alfvénic Mach number M A . The runs fall into two classes, that we will refer to throughout the paper, depending on the relative impact of turbulence and magnetic fields. Indeed, turbulence is subalfvénic in runs NoTurb (M A = 0) and SubA (M A < 1), and superalfvénic in runs SupA and SupAS (M A > 1). Regarding the Mach number, runs SupA and SubA have subsonic turbulence with M = 0.5 while run SupAS has a supersonic turbulence with 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), provided their mass at a given time and their accretion rate averaged over their lifetime. We model the dust sublimation by decreasing progressively the dust-to-gas ratio with the temperature. Gray opacities are taken from Semenov et al. (2003) and the dust sublimation is modelled by a dust-to-gas ratio which vanishes at high temperatures, similarly to Kuiper et al. (2010a). The dust-to-gas ratio varies as where M dust M gas 0 = 1% is the initial dust-to-gas mass ratio, and the evaporation temperature is where g = 2000 K and β = 0.0195 (Isella & Natta 2005). The gas opacity is taken equal to 0.01 cm 2 g −1 , so the total opacity tends towards 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 is the minimal optical depth allowed numerically (10 −4 ). 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 modelled 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).

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 2 . 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 3 (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 are 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. The radius of this sphere, so-called accretion radius is set to four times the finest resolution, i.e. 20 AU (40 AU for the LR runs). We use a threshold accretion scheme Table 1. Initial conditions of the four runs: name, Mach number, Alfvénic Mach number, mass-to-flux to critical mass-to-flux ratio, turbulence relative strength (with respect to thermal support and magnetic fields), respectively.

Model
M M A µ Turbulence relative strength NoTurb 0 0 5 No turbulence SupA 0.5 1.4 5 Superalfvénic, subsonic SupAS 2 5.7 5 Superalfvénic, supersonic SubA 0.5 0.57 2 Subalfvénic, subsonic based on the Jeans density in order to avoid artificial fragmentation. If a sink cell has a density above the local Jeans density, 10% of the excess is accreted by the sink particle within one time step. All of the accreted mass goes into 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 is prevented from diverging by the use of a so-called softening length, that is set equal to the accretion radius. Finally, sink particles can merge if their accretion radii overlap. We do not add any criterion to prevent sinks from merging.

Disk identification
Primary disk properties are presented in Sect. 3.3. These are computed from the cell-by-cell disk selection, which relies on the criteria presented in Joos et al. (2012) -It is rotationally-supported: ρv 2 φ /2 > f thres P, where v φ is the azimuthal velocity and P is the thermal pressure. We choose f thres = 2 as in Joos et al. (2012); -The gas number density is higher than n = 10 9 cm −3 , i.e. ρ = 3.85 × 10 −15 g cm −3 ; -The gas is not about to free-fall onto the central object too rapidly: v φ > f thres v r , where v r is the radial velocity; -The disk vertical structure is in hydrostatic equilibrium: v φ > f thres v z , where v z is the vertical velocity.
Let us note that there is no connectivity criterion, hence the extremity of a high-density spiral arm can be considered as part of the disk while the inter-arm low-density region (due to the gas being swept) may not be. We compute 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 add a geometrical criterion to avoid perturbations from companions and their possible disks: the cell must be located less than 0.9 times the binary separation. This is justified a posteriori as individual disks are observed around each star.

Results
We first summarize some of the main and common features or 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).

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 perpendicularly 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 we get a large filament-like structure due to the stronger inner turbulent support (see below), the dense region (ρ 10 −16 g cm −3 ) becomes rapidly concentrated 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 attributed to the toroidal magnetic pressure not being able to launch an outflow because of turbulence. Indeed, 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 and see that the central region in run SupAS (third column, second row) is magnetically-dominated (β < 1) on 2000 AU scales while thermally-dominated (β > 1) matter is infalling. This illustrates the importance to accurately account for the coupling between 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 subalfvénic, no secondary sink forms. With superalfvénic turbulence (runs SupA and SupAS), a secondary long-lived sink particle forms in the primary sink accretion disk. We will study in more details the stellar multiplicity in Sect. 3.2, and the secondary and circumbinary disks properties in Sect. 3.4 and Sect. 3.5, respectively. When not mentioned otherwise, we will refer to the primary sink and to the primary disk, for conciseness.

Magnetic field evolution
As collapse occurs, the magnetic field strength is expected to increase in the central regions. The density-magnetic field strength histograms for runs (from left to right) NoTurb, SupA, SupAS and SubA at t = 50 kyr are shown in Fig. 3. At densities below ∼10 −15 g cm −3 , we recover the ideal MHD limit where 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 plateaulike feature is present in the four runs, and contrasts with ideal MHD calculations, as shown in the low-mass ) and high-mass (C21) regimes. This is due to ambipolar diffusion, which becomes dominant above ρ 10 −15 g cm −3 . The diffusion coefficient varies non-linearly with the magnetic field strength η AD ∝ B 2 /ρ which explains its strong regulating effect. The plateau is located between ∼0.1 G in the superalfvénic runs (SupA, SupAS), and ∼0.3 G in the subalfvénic runs (NoTurb, SubA). The inclusion of ambipolar diffusion prevents the magnetic field strength to increase, which would change the disk structure and possibly the outflows since a strong magnetic field is reasonably expected in the magneto-centrifugal mechanism (C21). The large dispersion observed in the bottom-left panel of Fig. 3 at low density is provoked by turbulence in the core.

Asymmetries
While the numerical setup in run NoTurb is initially axisymmetric and symmetric with respect to the (x = 0, y, z)−plane (that we Fig. 1. Density slices perpendicular (first and second row, x10 zoom) and parallel (third and fourth row, x5 zoom) 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: run NoTurb, SupA, SupAS, SubA. A mass density ρ = 10 −19 g cm −3 corresponds to a particle density n = 2.6 × 10 4 cm −3 and ρ = 10 −11 g cm −3 to n = 2.6 × 10 12 cm −3 . White dots represent sink particles.
will refer to as the "north-south" symmetry), these symmetries are broken. First, we observe that 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. We find that the timescale associated with the interchange instability is indeed small enough (compared to the local advection timescale) to justify the interchange instability as a good candidate. This phenomenon is not the only asymmetry arising in the simulation. Indeed, we observe the presence of filamentary structures linking the densest regions (where the sink-disk system is) to the envelope, that we will refer to as "streamers". These 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 −15 g cm −3 . They appear as a path for the accretion flow and pull the magnetic field lines along with them, which in turn form an hour-glass shape. Since 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 magneto-centrifugal jets. These streamers form perpendicularly 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 stream- ers. In runs SupA and SubA, the streamers are much thicker. This could be a hint of turbulent diffusion but it 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 seem 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 outflows reported in Wu et al. (2004) are monopolar and our aim is to study a nonidealized case that could be compared to observations.

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 (similarly to Klassen et al. 2016) and exert gravitational A&A proofs: manuscript no. aanda  for the outflows (using the selection criteria presented in Paper II), for the gravitational torque and for the magnetic torque. The first and third integrals are performed over the surface S of a cylinder centered onto the primary sink, of radius R = 1000 AU and height H = 1000 AU and is 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, since we expect gravitational torques to transport angular momentum radially rather than vertically. In order to have comparable values from one run and one timestep to the other, we divide those fluxes by the mass within the cylinder. Figure 4 shows F out /M, F mag /M and F r grav /M as a function of the sink age. Angular momentum transported by outflows is slightly larger in run NoTurb than in SupA and SubA. We find larger 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 smaller in run 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. On the right panel of Fig. 4, we observe that the gravitational torque is stronger with increasing turbulence. Overall, the magnetic torque generalles dominates the angular momentum transport, except in run SupAS at later times where magnetic and gravitational torques are comparable. We perform the same comparison with R = 100 AU and the results remain qualitatively unchanged.

Sink mass history and stellar multiplicity
To begin with, we study the mass evolution of the sink particles. Left panel of Fig. 5 displays the primary and secondary (when there is) 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 behaviours are visible, between subalfvénic (NoTurb, SubA) and superalfvé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 subalfvénic cases than in the superalfvénic cases. This is confirmed when looking at 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, respectively. It is mainly between 10 −4 and 10 −3 M yr −1 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 most of the time zero values. We recall that our sink accretion scheme relies on the presence of high-enough density and Jeans-unstable gas within the sink volume. Hence, the absence of accretion, at a given time, means that the sink volume has not gathered enough material to be accreted. The main accretion events in the superalfvénic runs are more dramatic than in the subalfvénic runs, with companion sink particles or orbiting massive clumps raising the primary sink mass by a fraction of a solar mass instantaneously. Averaged over time, we obtain accretion rates in agreement with observational values  and references therein).
We report the formation of a long-lived binary system in the two superalfvénic runs (SupA and SupAS, see Table 2). In run SupAS, three additional sink particles form from initial fragmentation and four in the disk plane, but merge with the primary or secondary sinks. The secondary sink forms ≈ 17 kyr after the primary. It occurs at the extremity of a spiral arm.The secondary particle survives until the end of the run, i.e. a lifetime 33 kyr. This system also forms in the corresponding lowerresolution run (LRSupAS), at an age difference of ≈18 kyr instead of ≈17 kyr. Thanks to a lower resolution, run LRSupAS has been carried out up to t∼106 kyr and the secondary sink particle is ≈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 too 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 superalfvé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 the order of unity (in run LRSupA where it may tend towards 1), and even greater than 1 in run SupAS. It can be seen in the sink mass history (left panel of Fig. 5) that the primary sink is partially starved due to the presence of the secondary sink, as compared to the subalfvénic runs where 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 could be a more favorable place for accretion than the central location of the primary sink (when it dominates the total sink mass). Indeed, 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 formed from initial core fragmentation have merged. We recall that we use AMR based on the thermal Jeans length. The numerical convergence we perform with LR runs advocates for a physical fragmentation, rather than numerical. Nevertheless, the absence of a criterion to prevent sink merging in our simulations means that the final sink number may be higher. Hence, our multiplicity results can be considered, to first order, as lower-limit values. Fig. 6 shows the relation between the primary sink and the disk mass (left panel; the proper disk mass evolution is displayed in the right panel and 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 smaller in the subalfvénic runs. This is consistent with the gas falling smoothly onto the central star via the accretion disk, in subalfvénic runs, while in superalfvénic runs, clumps form in the disk and are subsequently accreted.

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 superalfvé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 most of the time between 50 and 200 AU in all runs. Large and ponctual increases coincide with the presence of a large spiral arm. The quasi-periodic variations found in runs SupA and SupAS are due to the orbital motions which impact the gas dynamics between the two stars.

Semi-analytical estimate of the disk radius
We compare the disk sizes with the theoretical predictions from Hennebelle et al. (2016a) for magnetically-regulated disk formation with ambipolar diffusion. Those are obtained from the equality between various timescales at the centrifugal radius. On the one hand, the timescale to generate toroidal field from differential rotation and the timescale for ambipolar diffusion to diffuse it vertically. On the other hand, the magnetic braking and the rotation timescales are set equal as well. The disk radius set where δ is the ratio between 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 further 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 (in opposite to the hydro case). Hence, the disk around primary sinks 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 analytically slightly overestimated for µ = 2 and µ = 5 runs, and underestimated when there is a significant rotational support (our SupA and SupAS runs, and the run MU5ADf in C21 with µ = 5 and E rot /E grav = 5%). This underestimation may be due to the presence of spiral patterns in disks when there is strong rotation/turbulent support, as the disk becomes gravitationally unstable.

Column density maps
Our disk definition is physically motivated (see Sect. 2.4) and allows for 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 be hardly checked observationally. For that purpose, Fig. 8 gives the column density maps at t = 50 kyr. The individual disks have a column density larger than their surroundings by at least one order of magnitude. Spiral arms are visible, as well as the presence of circumbinary gas. Figure 9 shows the azimuthal median of the column density as a function of the radius. Coloured circles indicate the radius which has been determined using our disk definition. We choose a median rather than a mean in order to reduce the impact of non-axisymetries due to fragments or sinks. A smooth but visible transition (in slope and values) is observed around the disk radius we have derived. For runs SupA and SupAS, the dense circumbinary gas makes more difficult the primary disk determination.

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 non-magnetic case studied in Mignon-Risse et al. (2020). Spiral arms are found in all runs, and at (nearly) all times. We check for the Toomre parameter (see Kratter & Lodato 2016 for a review on disk fragmentation) value before the development of spiral arms, as it indicates the unstable state to axisymmetric perturbations which lead to spiral arm development (often done in the literature, see e.g., Kratter & Matzner 2006, Klassen et al. 2016, Ahmadi et al. 2018. We compute it as in Mignon-Risse et al. (2020) -we do not display it here for conciseness. Disks are Toomre-unstable with typically Q 0.7, except in their outer parts (r 100 AU). Let us note that, as shown in Figs. 2 and 12, disks are (largely) dominated by thermal pressure rather than magnetic pressure (β > 1). Hence, computing the thermal Toomre (without the magnetic pressure) is sufficient.
However, the presence of spiral arms does not indicate fragmentation neither the formation of a multiple stellar system. In run NoTurb, no fragment forms. In run SubA, the first fragments form at a sink age around ∼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 before 16 kyr. These differences suggest that the interplay of turbulence and magnetic fields may impact disk fragmentation and sink formation.
Several criteria have been frequently used in the litterature 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 smaller 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. Therefore, cooling is not responsible for fragmentation here. Indeed, the Gammie model is well adapted to disks possessing local cooling processes while our disks heat and cool radiatively. Moreover, Gammie (2001) hypothesizes the disk to be cool and very thin, while we get 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 5 cm s −1 , r≈100 AU and Ω∼10 −11 s −1 , we get c s /(rΩ)∼10 0.12. Hence, we argue that the model of Gammie (2001) does not apply to the massive protostellar disks formed in this study.
We investigate spiral arm collision as a possibility to form fragments. Figure 10 shows the gas density as a function of position (along 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, i.e. a fragment, and decouples it from its parent arm. Collisions are observed in runs SubA and SupAS too. This process is favored when the spiral arms can sufficiently extend radially away from the primary star, which does not occur in run NoTurb. The rapid growth of the spiral arms in the turbulent runs can be seen in the left panel of Fig. 7. When turbulence is included, it brings additional angular momentum for spiral arms to extend (Hennebelle et al. 2016b, Hennebelle et al. 2017. The initial distribution of angular momentum computed with respect to the center of the domain increases with the distance. Thus, the less turbulent the core is, 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 disk fragmentation is due to spiral arm collision and favored by turbulence.
In order to understand the collapse of a fragment, let us follow the properties of the fragment of SupA presented above. Cells with a density larger than 10 −12 g cm −3 and a distance to the primary sink larger than 250 AU (to avoid the disk) are 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 sink formation time (t = 46.71 kyr). The free-fall time is computed as in Eq. 4 with the fragment mean density, and the sound-crossing time is equal to 2R/c s with 2R the estimate on the fragment diameter and c s its mean sound speed. This ratio gives an order-of-magnitude estimate of the fragment's stability to density perturbations: a ratio smaller than one indicates 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   We are particularly interested in the spiral arm collision occurring at t = 46.05 kyr, indicated by the two horizontal lines located at ±0.04 kyr from the collision time.
density are found roughly constant with time, the fragment radius increases by a factor ∼2.5 compared to its initial value. The fragment has accreted mass. Moreover, those quantities directly dictate the evolution of the ratio between the free-fall time and 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. sound-crossing time. Indeed, it is proportional to c s /( √ ρR), thus it tends towards an unstable state as R increases. Figure 11 shows 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 on fragment formation does not allow us to distinguish runs SubA and SupA, i.e. if there is an impact 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 a 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 can be explained by the magnetic braking being stronger in SubA than in SupA (see Fig.4 and Sect. 3.6). It takes longer for larger-scale gas -which carries more angular 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 subalfvenic, magnetic braking delays this process and drives the inward migration of fragments, preventing longlived stellar systems to emerge. 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 is in agreement with a Keplerian profile. It is slightly super-Keplerian in turbulent runs (SupA, SupAS, SubA) at radii 60 AU. In all runs, it becomes sub-Keplerian as the radius decreases. This is due to the gravitational field being dominated by the central object and diminished by the sink softening length mentioned in Sect. 2.3. The disks are roughly Keplerian, hence the infall velocity is much smaller than the free-fall velocity and typically smaller than 1 km s −1 . The rotation motions, and infall motions beyond the disk, are supersonic. The cells in the primary disk plane in the binary systems appear 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 a strong stellar activity from one component, they could be identified as large disks (Johnston et al. 2015).

Characteristic velocities and magnetic field components
As shown in the last row of Fig. 2, all primary disks have plasma beta β > 1(thermally-dominated). From density maps and Fig. 12 we observe that the disk radius is close to the point where a change from thermally-dominated to magneticallydominated region is observed. We argue that a physicallymotivated criterion for the identification of individual disks is the plasma beta with β 1 (equivalent to P th P mag ), since c s /v A = γβ/2∼ √ β, especially at late times. In our simulations, it encapsulates that the disk size is regulated by ambipolar diffusion, in contrast to the ideal MHD case (Masson et al. 2016, C21). This criterion only is not sufficient though because of the existence of thermally-dominated (β > 1) filaments, as well as 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 select cells in the disk plane but these are not restricted to the disk selection in order to probe the outer regions too. Strikingly, the vertical component B z dominates for r 50 AU in runs No-Turb and SubA, and in run SupA most of the computational time. A dominantly poloidal magnetic field is a necessary condition for launching 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 of B z and B r having opposite evolutions. Changes between B z and B r could be explained by the orbital motion of the primary sink, which occurs at superalfvenic 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 define the orbital Elsasser number for ambipolar diffusion Am as the ratio between the orbital time τ orb and the ion-neutral collision time t in = η AD c 2 /(4πv 2 A ) Taking the values in the vicinity of the primary sink we have v A ∼0.1km s −1 (see Fig. 12), η AD = 0.1, and τ orb 1 kyr (Fig. 7), which gives Am∼0.4 which is of the order of unity (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. This is due to the magnetic field lines being 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 which fan out at infinity and form an hourglass shape (e.g., Galli & Shu 1993).

Secondary disk (runs SupA, SupAS)
Let us focus 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 azimuthallyaveraged 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 where there is an inversion in the azimuthal velocity v φ corresponds to the closest part of the primary disk which 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, in complement to 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 (Sec. 3.3.5). The co-evolution of B z and B r likely results from the same mechanism, i.e. a B z component increased by ambipolar diffusion but turning into B r as the secondary sink moves at a superalfvenic speed. Nonetheless, the temporal evolution does not show a clear pattern permitting us to link the aforementionned observations 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, one may expect a magnetic structure favorable to outflow launching to build more rapidly around the secondary sink in run SupA than in run SupAS.
As can be seen on density maps (not shown here) and velocity profiles (e.g. taking the radius at which v φ drops as an upperlimit, in Fig. 14), the secondary disk radius is found to be of the order of ∼100 − 150 AU, in agreement with the transition radius β∼1. This is the same order of magnitude as found for the primary disks, and consistent within a factor of 2 with magnetic regulation (Eq. 10 predicts ∼200 AU in these cases).

Circumbinary disk (runs SupA, 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 investigate whether the binary separation is consistent with hydrodynamical disk radii, whose size is set by the centrifugal barrier from initial angular momentum conservation, rather than magnetic regulation. Nonetheless, this calculation will lead 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 rot = A&A proofs: manuscript no. aanda M c (R c Ω) 2 /2 = M c (J/M c ) 2 /(2R 2 c ), where R c is the its radius and J = r<R c r × ρv dV its angular momentum, centered onto the primary sink. Equalling the rotational energy and the gravitational energy 3GM c /(5R) (assuming a uniform density) we obtain the hydro disk radius The quantity J/M c is displayed in Fig. 16 as a function of the integration radius, here we take the value for R = R c . For run SupA, J/M c 8 × 10 21 cm 2 s −1 , hence r d,hy 300 AU, and J/M c 10 22 cm 2 s −1 so r d,hy 400 AU for run SupAS. These values roughly meet the binary separation (see Sec. 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 which 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 as the second disk grows. At t = 50 kyr, which is close to the birth epoch of the two secondary sinks, the circumbinary disk has mostly β ≥ 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 thermallysupported gas infall continues actively in run SupAS. The circumbinary disk does not appear as an isolated structure from the two individual disks, but the fate of this accreting system (in-dividual 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 No-Turb and SupA), as displayed in Fig. 13. Hence, the magnetic field geometry cannot be used here to discriminate between a unitary and a binary system.

What sets the 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. One may expect the disk normal 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 non-ideal MHD are accounted for. Since 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 fields orientation to tend towards a random distribution. Here we investigate 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 3 AU) reveals the disk orientation. Figure 16 shows the specific angular momentum j defined as j = J/M = 1 M r<R r×ρv dV as a function of the spatial scale R for three epochs: t = 0, 30, 50 kyr. We take r as radial vector with respect to the primary sink position, except for the first snapshot,  where there is no sink and we take the center of the box. We recall that each run has initially a tiny rotational support (1%) of solid-body rotation aligned with the 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 deg, not shown here for readability). Figure 16 shows the increase of the 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 4 AU) in runs SupA and SubA (left and right panels), but not in run SupAS where it is actually smaller. This means that the turbulent gas is in counter-rotation with respect to the initial solid-body rotation imposed. In the two runs with superalfvénic turbulence, SupA and SupAS (left and central panels of Fig. 16), the dominating components at disk scales (< 10 3 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.
Let us focus on the influence of magnetic fields. The left and right panels of Fig. 16 only differ by the magnetic field strength:  M r<R r × ρv dV as a function of 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 between the specific angular momentum j and the x-axis (the initial magnetic fields orientation) as a function of the sphere radius R for runs SupA (left), SupAS (middle) and SubA (right). The same timesteps 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) µ = 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, i.e. the transport of angular momentum outwards. It prevents disk formation perpendicularly to the magnetic field and favors configurations where 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 largescale magnetic field. The orientation of the angular momentum varies significantly with the scale considered and with time. We get similar results for sub-and supersonic turbulences as 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 velocity field changing 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 where the initial turbulence is subalfvénic, does not favor the alignment between j and B.
Overall, the disk normal in our simulations is misaligned (50 − 85 deg, 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 reached yet 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 it should first 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.

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 non-turbulent 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 non-ideal 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 nor the disk radius, up to M 15.8 M . Meanwhile, the influence of magnetic fields and ambipolar diffusion is undeniable. With highresolution simulations including ohmic dissipation (and no ambipolar diffusion), Kölligan & Kuiper (2018) obtain disks of up to 10 3 AU. While the Hall effect is expected to reduce or increase the disk size depending on the alignment between the rotation axis and the magnetic field axis, this suggests the strong regulating role of ambipolar diffusion (Hennebelle et al. 2016a) to be unique among non-ideal 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, which are reminiscent of the accretion channels found by Seifried et al. (2015) and the bridges in the study of (Kuffmeier et al. 2019, Kuffmeier et al. 2020). In agreement with Seifried et al. (2015), ram pressure dominates over magnetic pressure, but we further show that magnetic pressure is also dominated by thermal pressure. Those structures develop even in the non-turbulent run No-Turb, and they seem to be associated with magnetic fields effects rather than pure consequences from turbulence, as shown in Kuffmeier et al. (2019). We note though the the width of the accretion streamers seems to depend on the the initial turbulent level, the larger the turbulence, the wider the streamer. This effect could be a consequence of magnetic turbulent reconnection occurring in the envelop which increases with the turbulent level (e.g., Santos-Lima et al. 2013). The physical origin of these accretion streamer should be investigated in dedicated studies.
The present work confirms that disk formation does not occur preferentially perpendicularly to the core-scale magnetic fields, but its orientation is likely driven by turbulence (even subalfvénic). This is in agreement with the study of Machida et al. (2019) in the low-mass regime, in which they vary the angle between the rotation axis and the magnetic field direction. They observe 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. Such balanced mass ratios (of the order of unity) have been obtained in the radiationhydrodynamical 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 separation of 1590 AU. For comparison, binary separation is 350 − 600 AU in run SupA and 400−700 AU in run SupAS (the orbits are elliptic and the separation is slightly increasing 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 since they do only use a sink particle algorithm for the primary star.

Comparison with observations
Let us first compare our findings in terms of disk radius with observations. The disk radius of HH80-81, estimated to be ∼291 AU  or Orion Source 1 with 100 AU (Matthews et al. 2010), are in better agreement with the magnetically-regulated indivual disks 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 upperlimit of 125 AU set on the disk of S255IR SMA1 (Liu et al. 2020). Meanwhile, 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 massive protostars binary. Recent measurements with ALMA indicate the presence of jets, a circumbinary disk of radius ∼ 2500 AU and individual disks on scale ∼100 AU and 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 hints of counter-rotating disks, indicating preferentially another origin than disk fragmentation. Furthermore, Aizawa et al. (2020) have studied the disk orientation in five nearby star-forming regions and observe a random orientation, consistent with the disk plane being set by turbulence. Hence, it agrees with the present 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. 2019).

Limitations
Let us 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 have chosen a high-mass pre-stellar core for our initial conditions, consistent with the Turbulent Core model of McKee & Tan (2003). Indeed, this approach allows 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 leading to smaller disks than in the hydrody-