# Protoplanetary Disk Birth in Massive Star-forming Clumps: The Essential Role of the Magnetic Field

Ugo Lebreuilly, Patrick Hennebelle, Tine Colman, Benoît Commerçon, Ralf Klessen, Anaëlle Maury, Sergio Molinari, Leonardo Testi

► To cite this version:

Ugo Lebreuilly, Patrick Hennebelle, Tine Colman, Benoît Commerçon, Ralf Klessen, et al. Proto-planetary Disk Birth in Massive Star-forming Clumps: The Essential Role of the Magnetic Field. The Astrophysical Journal Letters, 2021, 917, 10.3847/2041-8213/ac158c. insu-03745701

HAL Id: insu-03745701

https://insu.hal.science/insu-03745701v1

Submitted on 1 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Protoplanetary Disk Birth in Massive Star-forming Clumps: The Essential Role of the Magnetic Field

Ugo Lebreuilly<fnref n="1" />, Patrick Hennebelle<fnref n="1" />, Tine Colman<fnref n="1" />, Benoît Commerçon<fnref n="2" />, Ralf Klessen<fnref n="3" />,<fnref n="4" />, Anaëlle Maury<fnref n="1" />, Sergio Molinari<fnref n="5" />, and Leonardo Testi<fnref n="1" />,<fnref n="6" />

Received 2021 May 28; revised 2021 July 14; accepted 2021 July 17; published 2021 August 12

## Abstract

Protoplanetary disks form through angular momentum conservation in collapsing dense cores. In this work, we perform the first simulations with a maximal resolution down to the astronomical unit (au) of protoplanetary disk formation, through the collapse of 1000 $M_\odot$ clumps, treating self-consistently both non-ideal magnetohydrodynamics with ambipolar diffusion as well as radiative transfer in the flux-limited diffusion approximation including stellar feedback. Using the adaptive mesh-refinement code RAMSES, we investigate the influence of the magnetic field on the disks properties with three models. We show that, without magnetic fields, a population dominated by large disks is formed that is not consistent with Class 0 disk properties as estimated from observations. The inclusion of magnetic field leads, through magnetic braking, to a very different evolution. When it is included, small <50 au disks represent about half the population. In addition, about 70% of the stars have no disk in this case, which suggests that our resolution is still insufficient to preserve the smaller disks. With ambipolar diffusion, the proportion of small disks is also prominent and we report a flat mass distribution around 0.01–0.1 $M_\odot$ and a typical disk-to-star mass ratios of ~$10^{-2}$–$10^{-1}$. This work shows that the magnetic field and its evolution plays a prominent role in setting the initial properties of disk populations.

Unified Astronomy Thesaurus concepts: Star formation (1569); Protoplanetary disks (1300); Magnetohydrodynamics (1964); Radiative transfer (1335); Interstellar medium (847); Molecular clouds (1072); Protostars (1302)

## 1. Introduction

Protoplanetary disks are a natural consequence of angular momentum conservation during the protostellar collapse. Planet formation in these disks not only depends on local quantities (density and temperature profile) but also on their global physical properties such as their total gas and dust masses or their size (see Testi et al. 2014, for a review on planet formation). Evolved disks, around Class II–III young-stellar objects (YSOs), are now observed in significant numbers at high resolution (see, for example, Andrews et al. 2018). They have typical sizes of 100 au (Andrews et al. 2009; Sanchis et al. 2021) but seem to be lacking the material to form giant planets (Manara et al. 2018), which suggests that those form early on. Unfortunately, the initial properties of protoplanetary disks are still poorly constrained. In contrast to Class II–III disks, Young Class 0–I disks are typically more compact and evolve over short lifetimes so they are difficult to observe. In addition, they are still deeply embedded in a massive envelope (André 2002) making their observation even more difficult. Following pioneering works, such as by Looney et al. (2000), who inferred the presence of a disk around Class 0 objects, recent surveys such as Continuum And Lines in Young ProtoStellar Objects (CALYPSO; Maury et al. 2019) and VLA/ALMA Nascent Disk and Multiplicity (VANDAM; Segura-Cox et al. 2018; Tobin et al. 2020) have started to probe them more extensively. Although the disk component is often unresolved or marginally resolved at these stages, Maury et al. (2019) have reported typical radii of about <50 au. Smaller disks could remain undetected because of insufficient resolution and studies such as by Yen et al. (2015) even argued that objects, such as the Class 0

protostar B335, could have a disk smaller than 10 au or even no disk at all.

On the theoretical side, our understanding of disk formation has gained in maturity over the two last decades. Historically, there has been an angular momentum problem in star formation. We know from observations that the specific angular momentum of dense cores is not conserved during star formation and must be redistributed by an efficient physical process. Magnetic field, one of the most promising candidates, has been extensively investigated as a possible solution, both using ideal (Allen et al. 2003; Hennebelle & Teyssier 2008; Joos et al. 2012) and non-ideal (Li et al. 2014; Machida et al. 2016; Masson et al. 2016; Wurster et al. 2016; Zhao et al. 2016; Vaytet et al. 2018; Hennebelle et al. 2020b; Wurster & Lewis 2020; Lee et al. 2021) magnetohydrodynamics (MHD). In both the ideal and non-ideal case, it was shown that magnetic braking would prevent the formation of large, massive, and unrealistic disks such as those observed in purely hydrodynamical simulations (see, for example, Machida & Matsumoto 2011). In fact, in some cases in ideal MHD it would totally prevent disk formation. Several studies (e.g., Dapp & Basu 2010; Tomida et al. 2015; Hennebelle et al. 2016) have shown that this so-called magnetic braking catastrophe could be solved by taking into account the role of diffusive processes, such as ambipolar diffusion, which reduces the braking efficiency at disk-like densities. In this case, a small disk with an initial size of ~20 au that would grow later on, is expected. It was also shown that other effects such as the magnetic field misalignment (Joos et al. 2012; Gray et al. 2018) or turbulence (Santos-Lima et al. 2012) could also reduce the efficiency of the magnetic braking and lead to the formation of small disks.

<footnote n="1">¹AIM, CEA, CNRS, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France; ugo.lebreuilly@cea.fr</footnote>

<footnote n="2">²Univ Lyon, Ens de Lyon, Univ Lyon1, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69007, Lyon, France</footnote>

<footnote n="3">³Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany</footnote>

<footnote n="4">⁴Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, INF 205, D-69120, Heidelberg, Germany</footnote>

<footnote n="5">⁵Istituto Nazionale di Astrofisica—IAPS, Via Fosso del Cavaliere 100, I-00133 Roma, Italy</footnote>

<footnote n="6">⁶ESO/European Southern Observatory, Karl-Schwarzschild-Strasse 2 D-85748 Garching bei München, Germany</footnote>
In past studies, most magnetized simulations have been performed to study star formation in the case of low-mass and isolated dense cores. In reality, these cores are connected to the large-scale environment inside molecular clouds and most stars (and disks) are not born in isolation, but rather within turbulent and magnetized complexes with masses of $10^2-10^3 M_\odot$ called clumps (see Urquhart et al. 2014; Elia et al. 2017, for studies of massive clumps at the galactic scale). A few studies have been dedicated to connect the giant molecular clouds to the disk scales by zooming in to a few individual stars/disks (Kuffmeier et al. 2017, 2019), but to statistically understand the formation of protoplanetary disks, and in the light of the recent progress in observing them in a population at very young stages, it is fundamental to systematically model the protostellar collapse starting from massive clumps all the way down to the disk scale. To date, only Bate (2018) investigated disk population formation in such clouds, although without including a magnetic field. Given the theoretical importance of the magnetic field in the isolated case, it is crucial to study its impact on a self-consistent disk population within these large scale environments.

In this work, we present a study of the disk populations resulting from the collapse of magnetized $1000 M_\odot$ clouds, in both the ideal and non-ideal MHD (with ambipolar diffusion) framework, taking into account the radiative stellar feedback. We present three collapse calculations with a maximum resolution of 1.2 au performed with the adaptive mesh-refinement (AMR) non-ideal MHD and radiative transfer code RAMSES. In Section 2, we describe our numerical methods. Then in Section 3 we report on the disk populations obtained from our three models. Finally, we summarize our main results in Section 4.

## 2. Methods

Our models are computed using the AMR (Berger & Oliger 1984) finite-volume code RAMSES (Teyssier 2002; Fromang et al. 2006) and its extension to radiative transfer in the flux-limited diffusion approximation (Commerçon et al. 2011, 2014), non-ideal MHD (Masson et al. 2012), and sink particles (Bleuler & Teyssier 2014).

Initially, we consider $1000 M_\odot$ uniform clumps of temperature 10 K and radius set according to the thermal-to-gravitational energy ratio $\alpha$

$$ \alpha = \frac{5 R_0 k_B T_0}{2 G M_0 \mu_g m_H}, \quad (1) $$

with $k_B$ being the Boltzmann constant, $G$ being the gravitational constant, $m_H$ being the hydrogen mass, and $\mu_g = 2.31$ being the mean molecular weight. We set $\alpha = 0.008$ and have an initial radius of $\sim 0.38$ pc and density of $3 \times 10^{-19} \text{ g cm}^{-3}$. This initial condition is fairly typical of the clumps observed in the Milky Way according to both the ATLASGAL (Urquhart et al. 2014) and HI-GAL (Elia et al. 2017) surveys. We note that we also set a gas adiabatic index $\gamma = 5/3$.

We set an initial turbulent velocity with a power-spectrum of $k^{-11/3}$, which corresponds to a Kolmogorov spectrum, and random phases. The rms of this velocity field is such that the initial Mach number is equal to 7, which corresponds to a turbulent-to-gravitational energy ratio of about 0.4.

In two models, we set a vertical and uniform initial magnetic field using the mass-over-flux to critical-mass-over-flux ratio

$\mu = 10$ such as

$$ \mu = \left( \frac{M_0}{\phi} \right) / \left( \frac{M}{\phi} \right)_c, \quad (2) $$

where $\left(\frac{M}{\phi}\right)_c = \frac{0.53}{\pi} \sqrt{5/G}$ (Mouschovias & Spitzer 1976). This corresponds to a magnetic field strength of $\sim 9.4 \times 10^{-5} \text{ G}$. In the non-ideal case, we only consider the effect of ambipolar diffusion, which is most probably the dominant non-ideal mechanism at the density range that we consider (Marchand et al. 2016), although uncertainties remains about the strength of the Hall effect. The value of the ambipolar resistivity is computed as a function of the temperature, density, and magnetic field intensity according to the table of Marchand et al. (2016). Similarly, the Planck and Rosseland opacities, used for the radiative transfer, are computed using the tables described by Vaytet et al. (2013).

To accurately follow the multiple scales of the clump, we use the AMR capability of RAMSES. The size of cell $\Delta x$ is given by $\Delta x = \frac{L_{\text{box}}}{2^\ell}$ where $\ell$ is the level of refinement. In our models, we consider a $\sim 1.53 \text{ pc box}$ and an initially uniform grid with $\ell_{\text{min}} = 7$ (corresponding to a 2460 au or 0.012 pc resolution) and we then refine the grid up to a maximum level $\ell_{\text{max}} = 18$ (1.2 au or $5.84 \times 10^{-6} \text{ pc}$) according to the Jeans length $\lambda_{\text{Jeans}}$ to impose

$$ \Delta x \le \frac{\lambda_{\text{Jeans}}}{N}, \quad (3) $$

where $N > 10$ to respect the Truelove et al. (1997) criterion and avoid artificial fragmentation.

In this work, we use sink particles (Bleuler & Teyssier 2014) to mimic the behavior of fully formed stars and avoid the numerical difficulties of increasing the numerical resolution enough to resolve them. We form sink particles when the density reaches $n_{\text{thre}} = 3 \times 10^{13} \text{ cm}^{-3}$ as in the standard case of Hennebelle et al. (2020b). Once a sink forms, it is placed at the position of the peak of the corresponding clump of density threshold $n_{\text{thre}}/10$ and the gas within a region of $4\Delta x$ (where $\Delta x$ is the cell size) is accreted if its density is above $n_{\text{thre}}/3$. After that, a fraction $C_{\text{acc}} = 0.1$ of the mass above $n_{\text{thre}}/3$ and within the accretion volume is attributed to the sink at each timestep. This corresponds to the fiducial value explored in Hennebelle et al. (2020b). As they show, the value of this parameter does impact the evolution of the disk and therefore should be explored in the future.

A star, i.e., a sink, of mass $M_{\text{sink}}$ and radius $R_*$ is a source of luminosity when accreting mass. This so-called accretion luminosity is expressed as

$$ L_{\text{acc}} = f_{\text{acc}} \frac{G M_{\text{sink}} \dot{M}_{\text{sink}}}{R_*}, \quad (4) $$

where $f_{\text{acc}} \le 1$ is the fraction of the accretion gravitational energy that is radiated away. In all our models, we consider $f_{\text{acc}} = 0.1$, which corresponds to the low value investigated by Hennebelle et al. (2020a). In their study, they also investigated other values but concluded that for a cloud similar to the ones considered here, the value of this parameter did not greatly affect the stellar mass spectrum. The star radius and its luminosity $L_{\text{int}}$ are computed using the models of Kuiper & Yorke (2013). Once the luminosities are computed, as in
<table><caption>**Table 1**
Summary of the Three Models</caption><thead><tr><th>Name</th><th>N<sub>star</sub></th><th>N<sub>isol</sub></th><th>N<sub>disks</sub></th><th>N<sub>syst</sub></th><th>t<sub>50M<sub>☉</sub></sub> (kyr)</th><th>t<sub>160M<sub>☉</sub></sub> (kyr)</th><th>M&#773;<sub>disk</sub> M<sub>☉</sub></th><th>R&#773;<sub>disk</sub> (au)</th></tr></thead><tbody><tr><td>IMHD</td><td>147</td><td>73</td><td>31</td><td>16</td><td>105</td><td>118</td><td>0.021</td><td>50.6</td></tr><tr><td>NIMHD</td><td>191</td><td>104</td><td>42</td><td>18</td><td>103</td><td>117</td><td>0.041</td><td>46.3</td></tr><tr><td>HYDRO</td><td>212</td><td>128</td><td>102</td><td>27</td><td>98</td><td>112</td><td>0.037</td><td>60</td></tr></tbody></table>

**Note.** From left to right: name of the model, number of stars (i.e., sinks), number of isolated/highly separated stars, number of disks and close multiple systems, $t_{50M_{☉}}$ and $t_{160M_{☉}}$, median disk mass and median disk radius. All the quantities with a bar are time averaged between $t_{50M_{☉}}$ and $t_{160M_{☉}}$, the other quantities are measured at $t_{160M_{☉}}$.

Hennebelle et al. (2020a), the corresponding energy ($L_{int} + L_{acc}$)dt is uniformly distributed to the sink cloud particles over the sink accretion volume (for more details on the cloud-in cell interpolation used, see Bleuler & Teyssier 2014).

### 3. Results

We now introduce our three models. First, two magnetized models were computed with (NIMHD) and without (IMHD) ambipolar diffusion. The third model (HYDRO), is similar but without any magnetic field. The main properties of the models are summarized in Table 1. We integrated them until the cumulative mass in sink particles is $160 M_{☉}$, which is 16% of the initial cloud mass. For the sake of comparison, averages are made between $t_{50M_{☉}}$ and $t_{160M_{☉}}$ where $50 M_{☉}$ and $160 M_{☉}$ have been accreted, respectively.

We define the number of isolated stars as the number of stars without any close neighbor. It implicitly includes stars that are in long-period multiple systems, the study of which is beyond the scope of this Letter. Our definition of a close system indeed only considers stars within less than 50 au (see the Appendix). We also define the number of systems as the number of primary stars that have at least one close neighbor.<fnref n="7" />

At the end of the calculation, NIMHD has 191 sinks/stars, 18 close multiple systems, 104 isolated sinks, and 42 disks are detected. HYDRO has 212 sinks, which includes 128 sinks in isolation, around which we find 102 disks. Finally, IMHD only has 147 sinks, including 73 in isolation, and we detect only 31 disks. The varying number of sinks in the three models suggests that fragmentation is significantly suppressed in IMHD but is also mostly recovered in NIMHD.

The mean disk masses of the three models are comparable although the value is slightly lower for IMHD for which it is ~0.021 $M_{☉}$, against ~0.041 $M_{☉}$ for NIMHD and 0.037 $M_{☉}$ for HYDRO. In terms of radii, HYDRO has the largest disks with a median radius of ~60 au. Both IMHD and NIMHD produce an extra population of small disks because of magnetic braking. For that reason, the median disk radius is smaller for the two magnetized models (~46.26 au for NIMHD and ~50.6 au for IMHD). In addition to reducing the typical size of the disks, the magnetic braking is in fact very damaging for the population of small disks. As previously mentioned, the IMHD model only has 31 disks at the end of the simulation, while NIMHD has 42. The ratio of disk-hosting stars is lower by a factor of almost two in the magnetized models. This suggests that, although more disks form through the course of the simulation, the magnetic braking strongly affects the very small disks (<10–15 au) to the point that some of them essentially become unresolved and disrupted at the grid scale. The smaller disks of

the magnetized runs indeed typically have their rotation axis aligned with the large-scale magnetic field, which maximizes the braking efficiency (see Joos et al. 2012; Gray et al. 2018, for studies of the impact of misalignment on disks properties). This critical effect should be investigated in future works with various initial field strengths. The values presented above are summarized in Table 1.

Figures 1, 2, and 3 show density slices for the NIMHD, IMHD, and HYDRO models, respectively. The disks are displayed at the end of the simulation, when $160 M_{☉}$, from the initial $1000 M_{☉}$, have been transformed into stars. In Figure 1, and more generally in the NIMHD model, the disks typically have small radii, among the six ones presented here (disks 1, 8, 40, 85, 121, and 140), half of them are smaller than 55 au and the others are smaller than 100 au. Their size is most likely controlled by magnetic self-regulation that favors the formation of small disks (Hennebelle et al. 2016). The disks of HYDRO and IMHD are on average larger, although for very different reasons. In the case of IMHD, this is only a slight difference and small disks are still often formed by the effect of the magnetic braking, but they are typically short-lived as no diffusive effect is counter-balancing the braking. Eventually, the population of large, hydro-like, disks survives better. For HYDRO, small disks can form, but they are statistically rarer as no efficient mechanism leads to angular momentum transport outside the disks. The 12 disks displayed in the respective panels of the two models (Figures 2 and 3) are again fairly typical and have radii of the order of 80–100 au (and up to ~200 au for the largest ones presented here).

In Figure 1, we added (bottom three panels) edge-on column density slices centered around sink 8 at various scales (0.58 pc, $3 \times 10^4$ au and 500 au). From the bottom-left panel, it appears quite clear that stars are not evenly spread in the cloud but rather concentrated within a cluster of about $1-2 \times 10^4$ au. As can be seen (top-middle panel and bottom-right panel), sink/disk 8 is not deeply embedded in this cluster and actually is quite isolated. Therefore, it is unsurprising that its disk, but also the ones around sink 121 and 140 (that are also quite isolated), appear to be mostly unperturbed and well organized. Conversely, some stars in the cluster (such as sink 40) are quite significantly perturbed by close interactions. Generally speaking, interactions between the disks and neighboring stars are significant for the three models. The example of sink 40 in the NIMHD model is striking, as it is part of the aforementioned cluster and highly perturbed by the other stars born in the same filamentary structure. Similar interactions and flybys are also happening in IMHD and HYDRO (e.g., disk 133 for IMHD or disk 13, 39 and 109 for HYDRO). Among the three models, HYDRO shows more of these interactions for mainly two reasons: (i) because disks are larger and therefore are affected

<footnote n="7">⁷ By definition, $N_{star} \neq N_{isol} + N_{syst}$.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** NIMHD model at $t_{160M_⊙} = 117$ kyr. Top six panels: density slices for various disks around primary stars. The colorbar and spatial extension of the six slices are the same for the sake of comparison and they are alternatively displayed edge-on or face-on. Finally, we display the sink number, the stellar and disk mass and the disk radius for each density snapshot. Bottom three panels: edge-on column density slices at various scales centered around sink 8. Circles represent sinks for all the snapshots and all of them are centered around the primary.</figcaption></figure>

by flybys over larger distances, and (ii) because the HYDRO disks are fragmenting more efficiently.

All three models show a wide diversity of sub-structures and commonly observed (around more evolved YSOs; Andrews 2020) disk features such as spirals (often associated with the presence of

a companion, e.g., disk 11 and 72 in IMHD), warps (e.g., disk 1 in IMHD), accretion streamers (e.g., disk 121 in NIMHD), and circum-multiplicity (e.g., disk 109 in HYDRO). We also note that disks tend to be either strongly perturbed or relatively quiet and isolated (e.g., disk 155 in HYDRO). In the former case this is either because
<figure><img src="image_4.png" /><figcaption>Figure 2. IMHD at $t_{160M_{\odot}} = 118$ kyr. Same as the top six snapshots of Figure 1.</figcaption></figure>

of star-disk interaction or because disks accrete high density material (e.g., disk 85 in NIMHD, which is also a circumbinary disk according to our definition, as are disk 45 in IMHD and disks 39 and 109 in HYDRO). We note that we do not see clear outflows around most stars in our models, which is likely a consequence of (i) the high level turbulence that reduces the coherence of the flow and (ii) the lack of resolution that precludes the launching of disk winds and jets.

In the top panels of Figure 4, we show the distributions, normalized to have a y-axis that peaks at 1 (to compare their shapes), of the disk masses (left), the disk-to-star mass ratio (middle) and the disk radius (right) for the three models (denoted by the three colors). We add a vertical line that represents the minimum solar mass nebular (MSMN; Hayashi 1981), the information in terms of Jupiter mass ($M_{Jup}$) and the Class 0 disks mass as estimated by Tychoniec et al. (2020) and converted into gas mass assuming dust-to-gas ratio of 1% in the disk mass histogram. In the disk radius histogram we add the Class 0 disks distributions of the CALYPSO (Maury et al. 2019, and references therein) and VANDAM (Segura-Cox et al. 2018; Tobin et al. 2020) surveys merged together (in pink), we also add three vertical dotted lines that represent the median disk radius of the models. We stress that the initial conditions used in this work are likely too compact compared

to the clouds observed by these two surveys. However, so far they are the only statistical samples available at the Class 0 stage and therefore it is worth to make this comparison keeping in mind the possible bias. Finally, we also recompute a disk radius distribution from the mass of our primaries using the analytical estimate $r_{\text{ana}}$ from Hennebelle et al. (2016) that writes as

$$r_{\text{ana}} = 18 \text{ au } \delta^{2/9} \left(\frac{\eta_{\text{AD}}}{0.1\text{s}}\right)^{2/9} \left(\frac{B_z}{0.1\text{G}}\right)^{-4/9} \left(\frac{M_{\text{disk}} + M_{\text{sink}}}{0.1M_{\odot}}\right)^{1/3}, \quad (5)$$

where $\delta$ is a coefficient of the order of a few, $\eta_{\text{AD}}$ is the ambipolar resistivity, $B_z$ is the vertical magnetic field in the disk and $M_{\text{disk}}$ and $M_{\text{sink}}$ are the disk and star mass, respectively. For simplicity we assumed $\delta = 1$ and $\eta_{\text{AD}} = 0$. 1 s (as in the aforementioned study) while the other quantities are volume averaged within the disk. As a complementary information, the three middle panels of the figure show the non-normalized distributions that helps to compare the three models directly. The bottom panels show the cumulative distribution function of the disk radius with three different normalization (left: normalized as a probability; middle:
<figure><img src="image_5.png" /><figcaption>**Figure 3.** HYDRO at $t_{\parallel 60 M_{\odot}} = 112$ kyr. Same as the top six snapshots of Figure 1.</figcaption></figure>

normalized by the ratio of disk-hosting star; right: normalized by the number of disks).

The three mass distributions are quite different, although their typical values are similar and range between $10^{-2}$ and $10^{-1} M_{\odot}$. As could be expected, HYDRO forms the most massive disks while IMHD forms the less massive ones despite forming the most massive stars. In both magnetic runs, some large disks still form when the magnetic field is misaligned and the braking efficiency is low. For IMHD, particularly, they are in fact even more stable than their hydrodynamical counterpart because their rotation generate a strong toroidal magnetic field that stabilize them against fragmentation (as it is the case, for example, for disk 72 in Figure 2). Contrary to the two other models, NIMHD has a flat disk mass distribution between $\sim 5 \times 10^{-2} M_{\odot}$ and $\sim 10^{-1} M_{\odot}$. Let us now concentrate on the disk-to-star mass ratio histogram. In HYDRO, it peaks around 0.05–0.2, which is very similar to what Bate (2018) has found in a previous study. IMHD forms less massive disks relatively to their parent star because stars are typically more massive as the level of fragmentation of this model is reduced. Finally, the distribution of NIMHD lies in between the ones of HYDRO and IMHD and peaks around $\sim 0.03–0.2$. In all the cases, the disk-to-star mass ratio histograms are quite peaked, which shows that the mass of the star and the mass of the disk are correlated. A

notable spread of about one order of magnitude, however, indicates that the stellar mass is not the only parameter that controls the disk mass.

We now focus on the disk radius histogram. We see that the NIMHD and IMHD models both form a population of small disks ($<50$ au), which represents about half of the disks. The HYDRO distribution of the radius peaks close to $\sim 60–70$ au. HYDRO presents the largest fraction of very large disks, with a prominent tail for the distributions between 100 au and 300–400 au, whereas the large disk population of both NIMHD and IMHD typically sharply decreases around 100–200 au. Although the shapes of the three radius histograms are quite similar, the three distributions are in fact very different as the three models do not form the same number of disks. This is shown clearly in the middle and bottom panels of Figure 4 where we display non-normalized histograms (of the disk mass, the disk-to-star mass ratio, and the disk radius) and cumulative distributions (of the disk radius) with three different normalization, respectively. Because about 55% of the stars are without a disk or with an unresolved disk component in NIMHD and IMHD, this means that only $\sim 20%–25\%$ of the stars have a disk larger than 50 au in these runs. The existence of a disk around these stars is uncertain and should be studied at higher resolution, which is not possible with our current numerical capabilities. However, it is strongly suggested by high resolution studies such as for example by
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top row: distributions (normalized so that their peak in the y-axis is 1) for the three models. Top-left panel: mass of the disks in solar mass (bottom x-axis) and Jupiter mass (top x-axis), minimum solar mass nebula (Hayashi 1981; dashed line), VLA Class 0 disk mass distribution (purple) inferred by Tychoniec et al. (2020). Top-middle panel: disk-to-primary-mass ratio. Top-right panel: disk radius distribution, medians radii of the models (dashed vertical lines, the same color as the models), analytical estimate (Hennebelle et al. 2016; black line), Class 0 disks (pink) from the merged CALYPSO (Maury et al. 2019), and VANDAM (Segura-Cox et al. 2018; Tobin et al. 2020) surveys. Middle row: the same as the top row but with no normalization of the histograms and without the observations. Bottom row: cumulative distribution function of disk radius for the three models. Bottom-left panel: normalized. Bottom-middle panel: normalized and multiplied by the ratio disk-hosting stars. Bottom-right panel: normalized and multiplied by the number of disks. All distribution are time averaged (between $t_{50M_\odot}$ and $t_{160M_\odot}$) distributions.</figcaption></figure>

Dapp & Basu (2010), Vaytet et al. (2018), that very small disks could form around some stars. We note that the population of small disks is slightly larger in NIMHD, which is consistent with the regulation of magnetic braking by ambipolar diffusion. Larger statistical samples would be useful to better access this

difference. Comparing the NIMHD probability density function and the one obtained with the analytical model, we find a reasonable agreement between the two distributions, although the analytical one presents a narrower range and has a less pronounced small disk population.
We superimposed the distributions of the radius of Class 0 disks extracted from the CALYPSO and VANDAM surveys merged together. Their histogram typically peaks around 50–60 au and agrees best with the ones from the two magnetized models. We note that our tentative comparison relies on both the observational and numerical definitions of a disk that both need to be questioned. It also relies on the underlying models that are used to estimate the disks properties from observations (for example assuming 1% of dust, which is not necessarily correct; see for e.g., Lebreuilly et al. 2020) and on detection bias, as very small disks might be unresolved around Class 0 protostars. In addition, our clouds have different conditions as the ones probed by these surveys (e.g., in the Perseus and Orion region), we rather model more compact clouds that have more affordable simulation cost. Future comparisons with real disks should be dedicated to generate synthetic observations of the models in order to ensure that the observational and numerical definition of a disk are in agreement.

## 4. Conclusion

In this work, we presented the first collapse calculations of massive 1000 $M_\odot$ clouds having au spatial resolution, which include a full treatment of the radiative transfer with stellar feedback as well as non-ideal MHD with ambipolar diffusion. Through three calculations, we explored the impact of the magnetic field on the disk populations self-consistently formed in the calculations. Our main findings and conclusions are as follows.

1. In the three models, we extract disk populations and infer their disk mass and radius distributions.

2. In the hydrodynamical case, and in accordance with the previous study of Bate (2018), we form massive and radially extended disks, often prone to fragmentation. We also note that, on average, about 70% of systems (single or multiple) have a resolved disk in this model.

3. In both the non-ideal and ideal MHD runs, we report a population of small disks. With ambipolar diffusion, the population is slightly prominent as magnetic braking is regulated by ambipolar diffusion. For both models, we note that the magnetic braking is particularly damaging for the population of small disks, leading to a low time-averaged ratio of disk-hosting stars about 44%–45% in both models. As mentioned in Section 3, very small disks might still form around the remaining stars if those were not dissipated at the grid scale, i.e., with a higher resolution. This stresses the need for development of new numerical methods to meet the challenge of computing more resolved models.

4. The disks formed in our calculations are initially massive enough to host solar-like exoplanetary systems. More than half of the disks are more massive than the MSMN in all our models. This assertion is strengthened by the fact that the limit has been downward revised in more recent theories (Desch 2007).

5. We produce a wide diversity of structures in our disks such as spirals, warps that are typical of older Class II–III objects. We also frequently observe stellar encounters and close multiplicity. Future dedicated studies should investigate whether they are related or whether the ones around older objects come from a different origin.

6. Fragmentation is significantly reduced in the ideal MHD case but mostly retrieved when including ambipolar diffusion.

We thank the referee for providing useful comments that helped us to significantly improve our manuscript. This work is funded by the ERC synergy ECOGAL (grant 855130, PIs: P. Hennebelle R. Klessen, S. Molinari, L. Testi) and was granted access to the HPC resources of CINES (OCCIGEN) under the allocation DARI A0090407023 made by GENCI. This work was concluded during the core2disk 2 program supported by Institut Pascal at Université Paris-Saclay with the support of the program “Investissements d’avenir” ANR-11-IDEX-0003-01. The figures were produced with matplotlib (https:// matplotlib.org/) using the OSYRIS library (https://pypi.org/project/ osyris/) developed by Neil Vaytet whom we thank. The data presented in this article will be made available on the Galactica database (https://galactica-simulations.eu/db/) developed by Damien Chapon, whom we also thank.

## Appendix

### A Disk Finder

The disk selection process is a challenging task. Disks are poorly defined and have arbitrary geometries. In addition, they can easily be confused with dense free-falling filamentary material and stars are often born in multiple systems. In this work, we propose a new way to isolate at best the disk material, which is largely inspired from Joos et al. (2012), who used the following criteria for a gas cell to be considered being included in the disk:

1. it should rotate faster than it falls radially $v_\phi > 2v_r$,

2. it should rotate faster than it falls vertically $v_\phi > 2v_z$,

3. it should not be thermally supported $\frac{1}{2}\rho v_\phi^2 > 2P_{\text{th}}$, where $P_{\text{th}}$ is the thermal pressure,

4. it should be composed of dense material $n > n_{\text{thre}} = 10^9 \text{ cm}^{-3}$, where n is the gas number density.

In our case, the disk selection operates as follows.

1. We place ourselves in the co-moving frame of the barycenter of stars within a distance $r_{\text{system}} = 50$ au from the analyzed sink. Note that we only analyze disks around primary stars (the most massive star within $r < r_{\text{system}} = 50$ au). We then get the disk rotation axis by computing the direction of the angular momentum in a pre-selected region of radius $R_{\text{mom}} = 15$ au.

2. We select the cells that verify the Joos et al. (2012) criterion in the disk frame.

3. We aim to exclude the streamers, i.e., the large-scale free-falling material. We first select the region within the cylindrical radius $R_{05}$ that encloses 95% of the mass. We also do the same in the vertical direction. Then we decompose the disk in $N_\theta = 50$ regions around the polar axis in the mid-plane (within $4\Delta x$) to compute the local radius. At this stage, we skip the analysis if there are no disk elements in more than two-thirds of the regions. Otherwise, we compute the disk radius $R_{\text{disk}}$ as the azimuthal median of the local radii.

4. We re-apply the Joos et al. (2012) criterion in the regions where $r < R_{\text{thre}} = \min(2R_{\text{disk}}, R_{05})$. The arbitrary pre-factor of 2 allows to preserve some of the non-axisymmetric features (e.g., spirals).
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Schematic view of the four disk finder steps: (i) change of frame and coordinates, (ii) pre-selection, (iii) estimation of the radius, and (iv) refined final selection.</figcaption></figure>

In Figure 5, we show a schematic simplified view of the four
disk finder steps.

ORCID iDs

Ugo Lebreuilly https://orcid.org/0000-0001-8060-1890

Patrick Hennebelle https://orcid.org/0000-0002-0472-7202

Tine Colman https://orcid.org/0000-0002-2636-4377

Benoît Commerçon https://orcid.org/0000-0003-2407-1025

Ralf Klessen https://orcid.org/0000-0002-0560-3172

Anaëlle Maury https://orcid.org/0000-0002-3801-8754

Sergio Molinari https://orcid.org/0000-0002-9826-7525

Leonardo Testi https://orcid.org/0000-0003-1859-3070

References

Allen, A., Shu, F. H., & Li, Z.-Y. 2003, ApJ, 599, 351

André, P. 2002, EAS, 3, 1

Andrews, S. M. 2020, ARA&A, 58, 483

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL, 869, L41

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502

Bate, M. R. 2018, MNRAS, 475, 5618

Berger, M. J., & Oliger, J. 1984, JCoPh, 53, 484

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

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

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

Dapp, W. B., & Basu, S. 2010, A&A, 521, L56

Desch, S. J. 2007, ApJ, 671, 878
Elia, D., Molinari, S., Schisano, E., et al. 2017, MNRAS, 471, 100

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

Gray, W. J., McKee, C. F., & Klein, R. I. 2018, MNRAS, 473, 2124

Hayashi, C. 1981, PThPS, 70, 35

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

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Chabrier, G. 2020a, ApJ, 904, 194

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020b, A&A, 635, A67

Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25

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

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

Kuffmeier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7

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

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

Lee, Y.-N., Charnoz, S., & Hennebelle, P. 2021, A&A, 648, A101

Li, Z. Y., Banerjee, R., Pudritz, R. E., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 173

Looney, L. W., Mundy, L. G., & Welch, W. J. 2000, ApJ, 529, 477

Machida, M. N., & Matsumoto, T. 2011, MNRAS, 413, 2767

Machida, M. N., Matsumoto, T., & Inutsuka, S.-i. 2016, MNRAS, 463, 4246

Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3

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

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

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

Maury, A. J., André, P., Testi, L., et al. 2019, A&A, 621, A76

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

Sanchis, E., Testi, L., Natta, A., et al. 2021, A&A, 649, A19

Santos-Lima, R., de Gouveia Dal Pino, E. M., & Lazarian, A. 2012, ApJ, 747, 21

Segura-Cox, D. M., Looney, L. W., Tobin, J. J., et al. 2018, ApJ, 866, 161

Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 339

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

Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130

Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117

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

Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19

Urquhart, J. S., Moore, T. J. T., Csengeri, T., et al. 2014, MNRAS, 443, 1555

Vaytet, N., Chabrier, G., Audit, E., et al. 2013, A&A, 557, A90

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

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

Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037

Yen, H.-W., Takakuwa, S., Koch, P. M., et al. 2015, ApJ, 812, 129

Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050