<figure><img src="image_1.png" /></figure>

HAL
open science

# Collapse of turbulent massive cores with ambipolar diffusion and hybrid radiative transfer

R. Mignon-Risse, M. González, B. Commerçon, Joakim A Rosdahl

► To cite this version:

R. Mignon-Risse, M. González, B. Commerçon, Joakim A Rosdahl. Collapse of turbulent massive cores with ambipolar diffusion and hybrid radiative transfer: I. Accretion and multiplicity. Astronomy and Astrophysics - A&A, 2021, 652, pp.A69. 10.1051/0004-6361/202140617. cea-03319227

HAL Id: cea-03319227

https://cea.hal.science/cea-03319227

Submitted on 11 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Collapse of turbulent massive cores with ambipolar diffusion and hybrid radiative transfer

## I. Accretion and multiplicity

R. Mignon-Risse<fnref n="1" />,<fnref n="2" />✉, M. González<fnref n="1" />, B. Commerçon<fnref n="3" />, and J. Rosdahl<fnref n="4" />

Received 19 February 2021 / Accepted 30 May 2021

### ABSTRACT

**Context.** Massive stars form in magnetized and turbulent environments and are often located in stellar clusters. The accretion and outflows mechanisms associated with forming massive stars and the origin of the stellar multiplicity of their system are poorly understood.

**Aims.** We study the effect of magnetic fields and turbulence on the accretion mechanism of massive protostars and their multiplicity. We also focus on disk formation as a prerequisite for outflow launching.

**Methods.** We present a series of four radiation-magnetohydrodynamical simulations of the collapse of a massive magnetized, turbulent core of 100 $M_\odot$ with the adaptive-mesh-refinement code RAMSES, including a hybrid radiative transfer method for stellar irradiation and ambipolar diffusion. We varied the Mach and Alfvénic Mach numbers to probe sub- and super-Alfvénic turbulence and sub- and supersonic turbulence regimes.

**Results.** Sub-Alfvénic turbulence leads to single stellar systems, and super-Alfvénic turbulence leads to binary formation from disk fragmentation following the collision of spiral arms, with mass ratios of 1.1–1.6 and a separation of several hundred AU that increases with initial turbulent support and with time. In these runs, infalling gas reaches the individual disks through 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 have similar properties. We obtain mass accretion rates of $\sim 10^{-4} M_\odot \text{ yr}^{-1}$ onto the protostars and 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 rather than by magnetic fields. Even without turbulence, axisymmetry and north-south symmetry with respect to the disk plane are broken by the interchange instability and thermally dominated streamers, respectively.

**Conclusions.** Small ($\le 300$ AU) massive protostellar disks such as those that are frequently observed today can so far only be reproduced 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 for distinguishing streamers and individual disks from their surroundings.

**Key words.** accretion, accretion disks – stars: formation – magnetohydrodynamics (MHD) – turbulence – stars: massive – methods: numerical

## 1. Introduction

Massive stars ($>8 M_\odot$) are luminous and are one of the main sources of feedback at parsec and galactic scales, especially when they explode as supernovae (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 ($\ge 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

<footnote n="1">¹ AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, 91191 Gif-sur-Yvette, France
e-mail: raphael.mignon-risse@apc.in2p3.fr</footnote>

<footnote n="2">² Université de Paris, CNRS, AstroParticule et Cosmologie, 75013, Paris</footnote>

<footnote n="3">³ Univ Lyon, Ens de Lyon, Univ Lyon1, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69007, Lyon, France</footnote>

<footnote n="4">⁴ Univ Lyon, Univ Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230, Saint-Genis-Laval, France</footnote>
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 hampered by 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 and occur 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_\odot$ (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 of multiple 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 increase the 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 & Garatti et al. 2017).

Despite the lack of systematic constraints, the presence of disks around young massive protostars ($L < 10^5 L_\odot$) 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/\phi)/(M/\phi)_{crit}$ is the ratio of mass-to-flux to critical mass-to-flux, where $\phi$ is the magnetic flux. It indicates whether magnetic fields can ($\mu < 1$) or cannot ($\mu > 1$) prevent collapse on their own. Nevertheless, a $\mu \ge 1$ still affects the gas dynamics. Several studies agree on supercritical values of $\mu = 1–4$ (Falgarone et al. 2008) or even $\mu \sim 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 $\mu = 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 flash-light 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 uni-dimensional 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{aligned}
\frac{\partial \rho}{\partial t} + \nabla \cdot [\rho \mathbf{u}] &= 0, \\
\frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot [\rho \mathbf{u} \otimes \mathbf{u} + P_{\text{I}}] &= -\lambda \nabla E_{\text{fld}} + \frac{\kappa_{\text{P},*} \rho}{c} \mathbf{F}_{\text{M1}} + \mathbf{F}_{\text{L}} - \rho \nabla \phi, \\
\frac{\partial E_{\text{T}}}{\partial t} + \nabla \cdot [\mathbf{u} (E_{\text{T}} + P + B^2/2)] \\
&\quad - (\mathbf{u} \cdot \mathbf{B})\mathbf{B} - \mathbf{E}_{\text{AD}} \times \mathbf{B}] = -\mathbb{P}_{\text{fld}} \nabla \cdot [\mathbf{u} + \kappa_{\text{P},*} \rho c E_{\text{M1}} - \lambda \mathbf{u} \nabla E_{\text{r}} \\
&\quad + \nabla \cdot \left(\frac{c\lambda}{\rho \kappa_{\text{R,fld}}} \nabla E_{\text{r}}\right) - \rho \mathbf{u} \nabla \phi, \\
\frac{\partial E_{\text{M1}}}{\partial t} + \nabla \cdot \mathbf{F}_{\text{M1}} &= -\kappa_{\text{P},*} \rho c E_{\text{M1}} + \dot{E}_{\text{M1}}^*, \\
\frac{\partial \mathbf{F}_{\text{M1}}}{\partial t} + c^2 \nabla \cdot \mathbb{P}_{\text{M1}} &= -\kappa_{\text{P},*} \rho c \mathbf{F}_{\text{M1}}, \\
\frac{\partial E_{\text{fld}}}{\partial t} - \nabla \cdot \left(\frac{c\lambda}{\rho \kappa_{\text{R,fld}}} \nabla E_{\text{fld}}\right) &= \kappa_{\text{P,fld}} \rho c (a_{\text{R}} T^4 - E_{\text{fld}}), \\
\frac{\partial \mathbf{B}}{\partial t} - \nabla \times [\mathbf{u} \times \mathbf{B} + \mathbf{E}_{\text{AD}}] &= 0, \\
\nabla \cdot \mathbf{B} &= 0, \\
\Delta \phi &= 4\pi G \rho.
\end{aligned}
\tag{1}
$$

Here, $\rho$ is the material density, $\mathbf{u}$ is the velocity, $P$ is the thermal pressure, $\lambda$ is the FLD flux limiter, $E_{\text{fld}}$ is the FLD radiative energy, $\kappa_{\text{P},*}$ is the Planck mean opacity at the stellar temperature, $\mathbf{F}_{\text{M1}}$ is the M1 radiative flux, $\mathbf{F}_{\text{L}} = (\nabla \times \mathbf{B}) \times \mathbf{B}$ is the Lorentz force, $\phi$ is the gravitational potential, $E_{\text{T}}$ is the total energy $E_{\text{T}} = \rho\epsilon + 1/2\rho u^2 + 1/2B^2 + E_{\text{fld}}$ ($\epsilon$ is the specific internal energy), $E_{\text{M1}}$ is the M1 radiative energy, $\mathbf{B}$ is the magnetic field, $\mathbf{E}_{\text{AD}}$ is the ambipolar electromotive force, $\mathbb{P}_{\text{fld}}$ is the FLD radiative pressure, $\kappa_{\text{P,fld}}$ is the Planck mean opacity in the FLD module, $\kappa_{\text{R,fld}}$ is the Rosseland mean opacity, $a_{\text{R}}$ is the radiation constant, $\mathbb{P}_{\text{M1}}$ is the M1 radiative pressure, and $\dot{E}_{\text{M1}}^*$ 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_{AD} = \frac{\eta_{AD}}{B^2} [(\nabla \times \mathbf{B}) \times \mathbf{B}] \times \mathbf{B}, \qquad (2)$$

where $\eta_{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 $\kappa_{P,*}\rho c E_{M1}$ couples the M1 and the FLD methods through the equation of the evolution of the internal energy,

$$C_v \frac{\partial T}{\partial t} = \kappa_{P,*} \rho c E_{M1} + \kappa_{P,\text{fld}} \rho c (E_{\text{fld}} - a_R T^4). \qquad (3)$$

We used the ideal gas relation for the internal specific energy $\rho\epsilon = C_v T$, 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_\odot$) of the density profile is then

$$\tau_{\text{ff}} = \sqrt{\frac{3\pi}{32G\bar{\rho}_0}} \approx 24 \text{ kyr}, \qquad (4)$$

where G is the gravitational constant, and $\bar{\rho}_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)_0/(M/\Phi)_{\text{crit}}$, where $(M/\Phi)_0 = M_c/(\pi B_0 R_c^2)$ and $(M/\Phi)_{\text{crit}} = 0.53/(3\pi) \sqrt{5/G}$ (Mouschovias & Spitzer 1976). Strong ($\mu = 2, B_0 = 170 \mu G$) and moderate ($\mu = 5, B_0 = 68 \mu G$) magnetic fields were considered. A drawback of this uniform distribution is that the mass-to-flux ratio decreases as $\mu \sim 1/R$ and is higher in the inner parts of the core, with $\mu \approx 50$ in the central plateau (for runs with $\mu = 5$). This corresponds to a weakly magnetized medium. We expect the central magnetic field strength to increase as $B \propto \rho^{2/3}$ as the core contracts and before ambipolar diffusion starts to dominate (at $\rho \sim 10^{-15} \text{g cm}^{-3}$), 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) \propto 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_{\text{rot}}/E_{\text{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 $M$ and Alfvénic Mach number $M_A$. The

<table><caption>**Table 1.** Initial conditions of the four runs.</caption><thead><tr><th>Model</th><th><i>M</i></th><th><i>M</i><sub>A</sub></th><th><i>μ</i></th><th>Turbulence relative strength</th></tr></thead><tbody><tr><td>NOTURB</td><td>0</td><td>0</td><td>5</td><td>No turbulence</td></tr><tr><td>SUPA</td><td>0.5</td><td>1.4</td><td>5</td><td>Superalfvénic, subsonic</td></tr><tr><td>SUPAS</td><td>2</td><td>5.7</td><td>5</td><td>Superalfvénic, supersonic</td></tr><tr><td>SUBA</td><td>0.5</td><td>0.57</td><td>2</td><td>Subalfvénic, subsonic</td></tr></tbody></table>

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 ($M_A = 0$) and SUBA ($M_A < 1$), and super-Alfvé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 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), 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} \arctan\left(\frac{T - T_{\text{evap}}(\rho)}{100}\right)\right), \qquad (5)$$

where $(M_{\text{dust}}/M_{\text{gas}})_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 cm}^{-3}} \right)^{\beta}, \qquad (6)$$

where $g = 2000$ K and $\beta = 0.0195$ (Isella & Natta 2005). The gas opacity was taken to be equal to $0.01 \text{ cm}^2 \text{ g}^{-1}$, so that the total opacity tends toward this value as the temperature increases beyond $T_{\text{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^{-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 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}} \sim M_c/R_c^2$. In comparison, the acceleration due to the nearest (0.6 pc) periodic core is $\sim M_c/(3R_c)^2 = a_{\text{grav}}/9$. A cell effective resolution is $0.5^\ell$ (in units of the box length), where $\ell$ is the AMR level of refinement. The coarse grid is $32^3$
(level $\ell = 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. The radius 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_{\phi}^2/2 > f_{\text{thres}} P$, where $v_{\phi}$ is the azimuthal velocity and $P$ is the thermal pressure. We chose $f_{\text{thres}} = 2$ as in Joos et al. (2012). The gas number density is higher than $n = 10^9 \text{cm}^{-3}$, i.e. $\rho = 3.85 \times 10^{-15} \text{g cm}^{-3}$. The gas is not about to free-fall onto the central object too rapidly: $v_{\phi} > f_{\text{thres}} v_r$, where $v_r$ is the radial velocity. The disk vertical structure is in hydrostatic equilibrium: $v_{\phi} > f_{\text{thres}} v_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 \approx 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 ($\rho \ge 10^{-16} \text{g cm}^{-3}$) rapidly concentrates in a sphere of diameter $\sim 2000$ AU (center panels, second row of Fig. 1). This is reminiscent of the structure described by Machida & Hosokawa (2020), and we attributed it to the inability of the toroidal magnetic pressure to launch an outflow because of turbulence. We show in Fig. 2 the plasma beta ($\beta = P_{\text{th}}/P_{\text{mag}}$, where $P_{\text{th}}$ and $P_{\text{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 ($\beta < 1$) at 2000 AU scales while thermally dominated ($\beta > 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 particle forms 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 $\sim 10^{-15} \text{g cm}^{-3}$, we recover the ideal MHD limit within which $B$ increases with $\rho$. In runs NOTURB and SUBA, the high-$B$, low-$\rho$ 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 $\rho \ge 10^{-15} \text{g cm}^{-3}$. The diffusion coefficient varies nonlinearly with the magnetic field strength $\eta_{\text{AD}} \propto B^2/\rho$, which explains its strong regulating effect. The plateau is located between $\sim 0.1 \text{ G}$ in the super-Alfvénic runs (SUPA, SUPAS) and at $\sim 0.3 \text{ 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 asymmetries are broken. First, pockets of magnetized plasma ($\beta < 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
<figure><img src="image_2.png" /><figcaption>**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 \approx 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 $\rho = 10^{-19}$ g cm$^{-3}$ corresponds to the particle density $n = 2.6 \times 10^4$ cm$^{-3}$ and $\rho = 10^{-11}$ g cm$^{-3}$ to $n = 2.6 \times 10^{12}$ cm$^{-3}$. White dots represent sink particles.</figcaption></figure>

(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. Filamentary structures 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 ($\beta > 1$), unlike the gas that surrounds them. They have a density $\rho \gtrsim 10^{-15}$ g cm$^{-3}$. 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 outflows reported in Wu et al. (2004) are monopolar, and
<figure><img src="image_3.png" /><figcaption>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 $\beta > 1$ as well.</figcaption></figure>

our aim is to study a nonidealized case that could be compared to observations.

### 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}} = \left| \int_S \rho r v_{\phi} \mathbf{v} \cdot \mathrm{d}\mathbf{S} \right| \quad (7)$$

for the outflows (using the selection criteria presented in Paper II),

$$F_{\text{grav}}^{\text{r}} = \left| \int_{S_{\text{rad}}} r \frac{g_{\phi}}{4\pi G} \mathbf{g} \cdot \mathrm{d}\mathbf{S} \right| \quad (8)$$

for the gravitational torque, and

$$F_{\text{mag}} = \left| \int_S r \frac{B_{\phi}}{4\pi} \mathbf{B} \cdot \mathrm{d}\mathbf{S} \right| \quad (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_{\text{rad}}$ only accounts for the surfaces of this cylinder that are perpendicular to the cylindrical radial vector because we expect
<figure><img src="image_4.png" /></figure>

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

<figure><img src="image_5.png" /></figure>

**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.

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_{\text{out}}/M$, $F_{\text{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 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. 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.

## 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 \approx 15.8 M_{\odot}$ 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 SUBA. 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^{-4}$ and $10^{-3} M_{\odot} \text{yr}^{-1}$ in run NOTURB. The accretion rate over the primary sink in run SUPA, which includes initial turbulence, is
<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 2.** Sink masses at the end of the simulations.</caption><thead><tr><th>Model</th><th>t<sub>end</sub></th><th>M<sub>★,end</sub></th><th>M<sub>2,end</sub></th></tr></thead><tbody><tr><td>NOTURB</td><td>80.4</td><td>15.8</td><td>—</td></tr><tr><td>SUPA</td><td>71.5</td><td>8.8</td><td>7.7</td></tr><tr><td>SUPAS</td><td>78.5</td><td>6.2</td><td>9.9</td></tr><tr><td>SUBA</td><td>66.2</td><td>11.0</td><td>—</td></tr></tbody></table>

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 system also 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 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 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 than 1 in run SUPAS. The sink mass history shows (left panel of Fig. 5) that the primary sink is partially starved due to the presence of 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.

### 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
<figure><img src="image_7.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>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.</figcaption></figure>

variations in the super-Alfvénic runs. We obtain disk masses ranging from ≈1–8 $M_\odot$ (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_\odot$) 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.

### 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_{d,AD} \approx 18 \text{ AU} \times \delta^{2/9} \left(\frac{\eta_{AD}}{0.1\text{s}}\right)^{2/9} \left(\frac{B_z}{0.1\text{G}}\right)^{-4/9} \left(\frac{M_d + M_*}{0.1 M_\odot}\right)^{1/3}, \quad (10)$$

where $\delta$ 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 $\delta = 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 $\eta_{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 (in contrast 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 $\mu = 2$ and $\mu = 5$, and it is underestimated when rotational support is significant (our SUPA and SUPAS runs, and run MU5ADf in C21 with $\mu = 5$ and
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Column density maps at $t \approx 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.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

$E_{\text{rot}}/E_{\text{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.

### 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-axisymmetries 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 \lesssim 0.7$, except in their outer parts ($r \gtrsim 100$ AU). As shown in Figs. 2 and 12, disks are (largely) dominated by thermal pressure rather than magnetic pressure ($\beta > 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\Omega) \ll 0.12$ (Eq. (26) of Gammie 2001), where $\Omega$ is the orbital frequency. With $c_s \approx 10^5$ cm s$^{-1}$, $r \approx 100$ AU and $\Omega \sim 10^{-11}$ s$^{-1}$, we obtain $c_s/(r\Omega) \sim 10 \gg 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
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

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^{-12}$ g cm$^{-3}$ 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-fall time 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_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

<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

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 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 sub-Alfvénic, magnetic braking delays this process and drives the inward migration of fragments, preventing long-lived stellar systems to emerge.

### 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_\theta)$ plane. Overall, the azimuthal velocity agrees with a Keplerian profile. It is slightly super-Keplerian in turbulent runs (SUPA, SUPAS, and SUBA) at radii $\ge$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$^{-1}$. 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.,
<figure><img src="image_13.png" /><figcaption>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 disk mass). 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.</figcaption></figure>

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 $\beta > 1$ (thermally dominated). In density maps and Fig. 12 we observe that the disk radius is close to the point at which 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 $\beta \ge 1$ (equivalent to $P_{\text{th}} \ge P_{\text{mag}}$) because $c_s/v_A = \sqrt{\gamma\beta/2} \sim \sqrt{\beta}$, especially at late 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 ($\beta > 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 \lesssim 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 $\tau_{\text{orb}}$ and the ion-neutral collision time $t_{\text{in}} = \eta_{\text{AD}}c^2/(4\pi v_A^2)$,

$$ Am = \frac{4\pi v_A^2 \tau_{\text{orb}}}{c^2 \eta_{\text{AD}}}. \qquad (11) $$

Taking the values in the vicinity of the primary sink, we have $v_A \sim 0.1 \text{ km s}^{-1}$ (see Fig. 12), $\eta_{\text{AD}} = 0.1$, and $\tau_{\text{orb}} \approx 1$ kyr (Fig. 7), which gives $Am \sim 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_\phi$ 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).

### 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_\phi$ is inverted corresponds to the closest part of the primary disk that dominates the azimuthal average (seen also around the primary
<figure><img src="image_14.png" /><figcaption>**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**).</figcaption></figure>

<caption>**Fig. 14.** Same as Fig. 14 for the secondary sinks in runs SUPA (**left**) and SUPAS (**right**) when their mass is 5 $M_\odot$.</caption>

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 $\beta \sim 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_\phi$ 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 $\beta \sim 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).

### 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_c(R_c\Omega)^2/2 = M_c(J/M_c)^2/(2R_c^2)$, where $R_c$ is the radius and $J = \int_{r<R_c} \mathbf{r} \times \rho\mathbf{v}\,dV$ its angular momentum, centered on the primary sink. Equalling the rotational energy and the gravitational energy $3GM_c/(5R)$ (assuming a uniform density), we
<figure><img src="image_15.png" /><figcaption>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.</figcaption></figure>

obtain the hydro disk radius

$$r_{\text{d,hy}} \approx 100 \text{ AU} \times \left( \frac{J/M_c}{5 \times 10^{21} \text{ cm}^2 \text{ s}^{-1}} \right)^2 \left( \frac{M_c}{100 M_\odot} \right)^{-1} . \quad (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 \approx 8 \times 10^{21} \text{ cm}^2 \text{ s}^{-1}$, hence $r_{\text{d,hy}} \approx 300$ AU, and $J/M_c \approx 10^{22} \text{ cm}^2 \text{ s}^{-1}$ so $r_{\text{d,hy}} \approx 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 $\beta \ge 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.

## 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} \int_{R}^{\infty} \mathbf{r} \times \rho \mathbf{v} \, dV$ 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 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 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: $\mu = 5$ in run SUPA (left) and $\mu = 2$ in run SUBA (right). At small scales, the component $j_x$ in run SUPA is $\approx 2$ times larger than in
<figure><img src="image_16.png" /><figcaption>**Fig. 16.** Specific angular momentum **j** = **J**/*M* = (1/M) ∫<sub>Jr<R</sub> **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.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**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).</figcaption></figure>

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 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*<sub>y</sub> < *j*<sub>z</sub> (to be distinguished from *j*<sub>x</sub>, which is more affected by magnetic braking) at the disk scale, while *j*<sub>y</sub> > *j*<sub>z</sub> 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

($\le$ 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.

# 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*<sub>☉</sub>. This value is very similar to what has been found with the FLD method in C21 (their nonideal MHD run with *μ* = 5 and *E*<sub>rot</sub>/*E*<sub>grav</sub> = 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*<sub>☉</sub>. 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<sup>3</sup> 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 fur-
ther 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, how-
ever: the stronger the turbulence, the wider the streamer. This
effect could be a consequence of magnetic turbulent reconnec-
tion 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 ($\mu = 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_\odot$ and 29.2 $M_\odot$ and a separation of 1590 AU. For com-
parison, the binary separation is 350–600 AU in run SUPA and
400–700 AU in run SUPAS (the orbits are elliptic, and the sep-
aration 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 compa-
rable 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 obser-
vations. 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 magnet-
ically regulated individual 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, con-
sistent 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 hydrodynamami-
cal 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 gravi-
tational 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 struc-
ture, and conclude on the effect of protostellar heating on disk
fragmentation. A cell of gas density $10^{-11}$ g cm$^{-3}$ and opacity to
ultraviolet radiation ~$10^2$ cm$^2$ g$^{-1}$ would require a spatial resolu-
tion of less than $10^{-3}$ 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 vol-
ume 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 mix-
ture and radiation. Furthermore, dust grains are charge carriers,
which means that the ionization degree would vary and the non-
ideal 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 $\le 10^4$ 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.

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.

## References

Ahmadi, A., Beuther, H., Mottram, J. C., et al. 2018, A&A, 618, A46

Aizawa, M., Suto, Y., Oya, Y., Ikeda, S., & Nakazato, T. 2020, ApJ, 899, 55

Beltrán, M. 2020, ArXiv e-prints [arXiv:2005.06912]

Beltrán, M. T., & de Wit W. J. 2016, A&ARV, 24, 6

Beltrán, M. T., Padovani, M., Girart, J. M., et al. 2019, A&A, 630, A54

Beuther, H., Ahmadi, A., Mottram, J. C., et al. 2019, A&A, 621, A122

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015

Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 2001, MNRAS, 323, 785

Bosco, F., Beuther, H., Ahmadi, A., et al. 2019, A&A, 629, A10

Caratti o Garatti, A., Stecklum, B., Garcia Lopez, R., et al. 2017, Nat. Phys., 13, 276

Cesaroni, R., Neri, R., Olmi, L., et al. 2005, A&A, 434, 1039

Chandrasekhar, S., & Fermi, E. 1953, ApJ, 118, 113

Commerçon, B., Hennebelle, P., & Henning, T. 2011a, ApJ, 742, L9

Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011b, A&A, 529, A35

Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11

Commerçon, B., González, M., Mignon-Risse, R., Hennebelle, P., & Vaytet, N. 2021, A&A, submitted

Duchêne, G., & Kraus, A. 2013, Annu. Rev. Astron. Astrophys., 51, 269

Falgarone, E., Troland, T. H., Crutcher, R. M., & Paubert, G. 2008, A&A, 487, 247

Fernández-López, M., Girart, J. M., Curiel, S., et al. 2011, ApJ, 142, 97

Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43

Foglizzo, T., Scheck, L., & Janka, H.-T. 2006, ApJ, 652, 1436

Fontani, F., Commerçon, B., Giannetti, A., et al. 2016, A&A, 593, L14

Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371

Fuksman, J. D. M., Klahr, H., Flock, M., & Mignone, A. 2021, ApJ, 906, 78

Galli, D., & Shu, F. H. 1993, ApJ, 417, 220

Gammie, C. F. 2001, ApJ, 553, 174

Gaudel, M., Maury, A. J., Belloche, A., et al. 2020, A&A, 637, A92

Ginsburg, A., Bally, J., Goddi, C., Plambeck, R., & Wright, M. 2018, ApJ, 860, 119

Girart, J. M., Beltran, M. T., Zhang, Q., Rao, R., & Estalella, R. 2009, Science, 324, 1408

Girart, J. M., Frau, P., Zhang, Q., et al. 2013, ApJ, 772, 69

Girart, J. M., Fernandez-López, M., Li, Z.-Y., et al. 2018, ApJ, 856, L27

González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429

Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126

Hennebelle, P., & Fromang, S. 2008, A&A, 477, 9

Hennebelle, P., Commerçon, B., Joos, M., et al. 2011, A&A, 528, A72

Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016a, ApJ, 830, L8

Hennebelle, P., Lesur, G., & Fromang, S. 2016b, A&A, 590, A22

Hennebelle, P., Lesur, G., & Fromang, S. 2017, A&A, 599, A86

Hull, C. L. H., Gouvellec, V. J. M. L., Girart, J. M., Tobin, J. J., & Bourke, T. L. 2020, ApJ, 892, 152

Isella, A., & Natta, A. 2005, A&A, 438, 899

Johnston, K. G., Robitaille, T. P., Beuther, H., et al. 2015, ApJ, 813, L19

Joos, M., Hennebelle, P., & Giardi, A. 2012, A&A, 543, A128

Joos, M., Hennebelle, P., Giardi, A., & Fromang, S. 2013, A&A, 554, A17

Kaisig, M., Tajima, T., & Lovelace, R. V. E. 1992, ApJ, 386, 83

Keto, E. 2007, MNRAS, 666, 976

Klassen, M., Pudritz, R. E., Kuiper, R., Peters, T., & Banerjee, R. 2016, ApJ, 823, 28

Kölligan, A., & Kuiper, R. 2018, A&A, 620, A182

Krasnopolsky, R., Li, Z.-Y., Shang, H., & Zhao, B. 2012, ApJ, 757, 77

Kratter, K. M., & Matzner, C. D. 2006, MNRAS, 373, 1563

Kratter, K., & Lodato, G. 2016, Annu. Rev. Astron. Astrophys., 54, 271

Kraus, S., Hofmann, K.-H., Menten, K. M., et al. 2010, Nature, 466, 339

Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science, 323, 754

Kuffmeier, M., Calcutt, H., & Kristensen, L. E. 2019, A&A, 628, A112

Kuffmeier, M., Reissl, S., Wolf, S., Stephens, I., & Calcutt, H. 2020, A&A, 639, A137

Kuiper, R., & Yorke, H. W. 2013, ApJ, 772, 61

Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010a, ApJ, 722, 1556

*Acknowledgements.* This work was supported by the CNRS "Programme National de Physique Stellaire" (PNPS). The numerical simulations presented
R. Mignon-Risse et al.: Collapse of turbulent massive cores with ambipolar diffusion and hybrid radiative transfer. I.

Kuiper, R., Klahr, H., Dullemond, C., Kley, W., & Henning, T. 2010b, A&A, 511, A81

Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20

Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2012, A&A, 537, A122

Lam, K. H., Li, Z.-Y., Chen, C.-Y., Tomida, K., & Zhao, B. 2019, MNRAS, 489, 5326

Larson, R. B. 1974, MNRAS, 169, 229

Larson, R. B., & Starrfield, S. 1971, A&A, 13, 190

Lebreuilly, U., Commerçon, B., & Laibe, G. 2019, A&A, 626, A96

Lebreuilly, U., Commerçon, B., & Laibe, G. 2020, A&A, 641, A112

Levermore, C. D. 1984, JQSRT, 31, 149

Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321

Li, H.-b., Yuen, K. H., Otto, F., et al. 2015, Nature, 520, 518

Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180

Liu, S.-Y., Su, Y.-N., Zinchenko, I., et al. 2020, ApJ, 904, 181

Lovelace, R. V. E., & Scott, H. A. 1981, Plasma Astrophysics, Course and Workshop, 215

Machida, M. N., & Hosokawa, T. 2020, MNRAS, 499, 4490

Machida, M. N., Hirano, S., & Kitta, H. 2019, MNRAS

Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18

Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24

Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016, A&A, 587, A32

Matthews, L. D., Greenhill, L. J., Goddi, C., et al. 2010, ApJ, 708, 80

Maury, A. J., Girart, J. M., Zhang, Q., et al. 2018, MNRAS

McKee, C. F., & Tan, J. C. 2003, ApJ, 585, 850

Meyer, D. M.-A., Kuiper, R., Kley, W., Johnston, K. G., & Vorobyov, E. 2018, MNRAS, 473, 3615

Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2020, A&A, 635, A42

Motte, F., Bontemps, S., & Louvet, F. 2018, Annu. Rev. Astron. Astrophys., 56, 41

Mouschovias, T. C., & Spitzer, Jr. L. 1976, ApJ, 210, 326

Myers, A. T., McKee, C. F., Cunningham, A. J., Klein, R. I., & Krumholz, M. R. 2013, ApJ, 766, 97

Nony, T., Louvet, F., Motte, F., et al. 2018, A&A, 618, L5

Oliva, G. A., & Kuiper, R. 2020, A&A, 644, A41

Padoan, P., Pan, L., Juvela, M., Haugbølle, T., & Nordlund, Å. 2020, ApJ, 900, 82

Palau, A., Fuente, A., Girart, J. M., et al. 2013, ApJ, 762, 120

Patel, N. A., Curiel, S., Sridharan, T. K., et al. 2005, Nature, 437, 109

Peters, T., Banerjee, R., Klessen, R. S., & Low, M.-M. M. 2011, ApJ, 729, 72

Pillai, T., Kauffmann, J., Tan, J. C., et al. 2015, ApJ, 799, 74

Pillai, T., Kauffmann, J., Wiesemeyer, H., & Menten, K. M. 2016, A&A, 591, A19

Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nat. Astron., 4, 1158

Ramsey, J. P., & Dullemond, C. P. 2015, A&A, 574, A81

Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188

Rosdahl, J., & Teyssier, R. 2015, MNRAS, 449, 4380

Rosen, A. L., Krumholz, M. R., McKee, C. F., & Klein, R. I. 2016, MNRAS, 463, 2553

Rosen, A., Krumholz, M., Oishi, J., Lee, A., & Klein, R. 2017, J. Comput. Phys., 330, 924

Rosen, A. L., Li, P. S., Zhang, Q., & Burkhart, B. 2019, ApJ, 887, 108

Santos-Lima, R., de Gouveia Dal Pino, E. M., & Lazarian, A. 2013, MNRAS, 429, 3371

Schneider, N., Csengeri, T., Bontemps, S., et al. 2010, A&A, 520, A49

Seifried, D., Banerjee, R., Klessen, R. S., Duffin, D., & Pudritz, R. E. 2011, MNRAS, 417, 1054

Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2012, MNRAS Lett., 423, L40

Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2015, MNRAS, 446, 2776

Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611

Shu, F. H. 1977, ApJ, 214, 488

Tan, J. C., Beltran, M. T., Caselli, P., et al. 2014 Protostars and Planets VI, eds. H. Beuther, R. Klessen, C. Dullemond, Th. Henning (University of Arizona Press) 149

Tanaka, K. E. I., Zhang, Y., Hirota, T., et al. 2020, ApJ, 900, L2

Tang, Y.-W., Ho, P. T. P., Girart, J. M., et al. 2009, ApJ, 695, 1399

Teyssier, R. 2002, A&A, 385, 337

Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179

Vaytet, N., Commerçon, B., Masson, J., González, M., & Chabrier, G. 2018, A&A, 615, A5

Vázquez-Semadeni, E., Gómez, G. C., Jappsen, A.-K., Ballesteros-Paredes, J., & Klessen, R. S. 2009, ApJ, 707, 1023

Vázquez-Semadeni, E., González-Samaniego, A., & Colín, P. 2016, MNRAS, stw3229

Wu, Y., Wei, Y., Zhao, M., et al. 2004, A&A, 426, 503

Wurster, J., & Lewis, B. T. 2020, MNRAS, 495, 3807

Yorke, H. W., & Sonnhalter, C. 2002, ApJ, 569, 846

Zhang, Q., Qiu, K., Girart, J. M., et al. 2014, ApJ, 792, 116
## Appendix A: Luminosity injection into the sink particle volume: disk size

<figure><img src="image_18.png" /><figcaption>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.</figcaption></figure>

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).

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{\left\| \frac{\partial((\gamma - 1)P + B_s^2/2)}{\partial y} \right\|}{\left( (\gamma - 1)c_s^2 + v_A^2 \right) \left\| \frac{\partial \rho}{\partial y} \right\|} > 1, \qquad (B.1) $$

where x is the normal direction to the disk, and $v_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 $\omega = \text{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}}, \qquad (B.2) $$

where we defined that $\alpha \equiv v_A^2/c_s^2$, $s = \frac{1}{\gamma-1} \ln(P\rho^{-\gamma})$ is the normalized gas entropy, and $g_{\text{eff}} = g - v_{\phi}^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^2$ 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 $\omega \approx 70 \text{ kyr}^{-1}$ 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 $\tau_{\text{instab}} \approx 14 \text{ yr}$. The unstable zone is $\approx 20 \text{ AU}$ wide, in which the gas is flowing at a radial velocity $v_r \approx 0.8 \text{ km s}^{-1}$, so that the advection timescale is $\tau_{\text{adv}} \approx 120 \text{ yr}$. Because $\tau_{\text{instab}} \le \tau_{\text{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 though the interchange instability is a good candidate.
<figure><img src="image_19.png" /><figcaption>**Fig. B.1.** Interchange instability in run NOTURB. Top left and right panels: plasma $\beta$ at $t = 45.25$ kyr (before the instability develops) and $t = 45.37$ kyr. Bottom panel: square root $\omega$ of $N^2$ (see Eq. (B.2)): $\text{Im}(\omega)$ 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.</figcaption></figure>