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

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.

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²–10³ M_⊙ 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_⊙ 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.

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_⊙ uniform clumps of temperature 10 K and radius set according to the thermal-to-gravitational energy ratio α

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{5}{2}\displaystyle \frac{{R}_{0}{k}_{{\rm{B}}}{T}_{0}}{{ \mathcal G }{M}_{0}{\mu }_{{\rm{g}}}{m}_{{\rm{H}}}},\end{eqnarray} \tag{ 1 }$$

with k_B being the Boltzmann constant, ${ \mathcal G }$ being the gravitational constant, m_H being the hydrogen mass, and μ_g = 2.31 being the mean molecular weight. We set α = 0.008 and have an initial radius of ∼0.38 pc and density of 3 × 10⁻¹⁹ g cm⁻³. 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 γ = 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 μ = 10 such as

$$\begin{eqnarray}&&\mu =\left(\displaystyle \frac{{M}_{0}}{\phi }\right)/{\left(\displaystyle \frac{M}{\phi }\right)}_{c},\end{eqnarray} \tag{ 2 }$$

where ${\left(\tfrac{M}{\phi }\right)}_{c}=\tfrac{0.53}{\pi }\sqrt{5/{ \mathcal G }}$ (Mouschovias & Spitzer 1976). This corresponds to a magnetic field strength of ∼9.4 × 10⁻⁵ 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 Δx is given by ${\rm{\Delta }}x=\tfrac{{L}_{\mathrm{box}}}{{2}^{{\ell }}}$ where ℓ is the level of refinement. In our models, we consider a ∼1.53 pc box and an initially uniform grid with ${{\ell }}_{\min }=7$ (corresponding to a 2460 au or 0.012 pc resolution) and we then refine the grid up to a maximum level ${{\ell }}_{\max }=18$ (1.2 au or 5.84 × 10⁻⁶ pc) according to the Jeans length λ_Jeans to impose

$$\begin{eqnarray}&&{\rm{\Delta }}x\leqslant \displaystyle \frac{{\lambda }_{\mathrm{Jeans}}}{N},\end{eqnarray} \tag{ 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_thre = 3 × 10 ¹³ cm⁻³ 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_thre/10 and the gas within a region of 4Δx (where Δx is the cell size) is accreted if its density is above n_thre/3. After that, a fraction C_acc = 0.1 of the mass above n_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_sink and radius R_⋆ is a source of luminosity when accreting mass. This so-called accretion luminosity is expressed as

$$\begin{eqnarray}&&{L}_{\mathrm{acc}}={f}_{\mathrm{acc}}\displaystyle \frac{{ \mathcal G }{M}_{\mathrm{sink}}{\dot{M}}_{\mathrm{sink}}}{{R}_{\star }},\end{eqnarray} \tag{ 4 }$$

where f_acc ≤ 1 is the fraction of the accretion gravitational energy that is radiated away. In all our models, we consider f_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_int are computed using the models of Kuiper & Yorke (2013). Once the luminosities are computed, as in 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).

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{| }_{50{M}_{\odot }}$ and $t{| }_{160{M}_{\odot }}$ where 50M_⊙ and 160M_⊙ have been accreted, respectively.

Table 1. Summary of the Three Models

Name	Nstar	Nisol	Ndisks	Nsyst	$t{| }_{50{M}_{\odot }}$ (kyr)	$t{| }_{160{M}_{\odot }}$ (kyr)	${\overline{M}}_{\mathrm{disk}}$ M⊙	${\overline{R}}_{\mathrm{disk}}$ (au)
IMHD	147	73	31	16	105	118	0.021	50.6
NIMHD	191	104	42	18	103	117	0.041	46.3
HYDRO	212	128	102	27	98	112	0.037	60

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{| }_{50{M}_{\odot }}$ and $t{| }_{160{M}_{\odot }}$, median disk mass and median disk radius. All the quantities with a bar are time averaged between $t{| }_{50{M}_{\odot }}$ and $t{| }_{160{M}_{\odot }}$, the other quantities are measured at $t{| }_{160{M}_{\odot }}$.

Download table as:  ASCIITypeset image

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

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 1, we added (bottom three panels) edge-on column density slices centered around sink 8 at various scales (0.58 pc, 3 × 10⁴ 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 × 10⁴ 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 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 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_ana from Hennebelle et al. (2016) that writes as

$$\begin{eqnarray}\begin{array}{rcl}{r}_{\mathrm{ana}} & = & 18\,\mathrm{au}\,{\delta }^{2/9}{\left(\displaystyle \frac{{\eta }_{\mathrm{AD}}}{0.1{\rm{s}}}\right)}^{2/9}{\left(\displaystyle \frac{{B}_{z}}{0.1{\rm{G}}}\right)}^{-4/9}\\ & & {\left(\displaystyle \frac{{M}_{\mathrm{disk}}+{M}_{\mathrm{sink}}}{0.1{M}_{\odot }}\right)}^{1/3},\end{array}\end{eqnarray} \tag{ 5 }$$

where δ is a coefficient of the order of a few, η_AD is the ambipolar resistivity, B_z is the vertical magnetic field in the disk and M_disk and M_sink are the disk and star mass, respectively. For simplicity we assumed δ = 1 and η_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: normalized by the ratio of disk-hosting star; right: normalized by the number of disks).

Zoom In Zoom Out Reset image size

The three mass distributions are quite different, although their typical values are similar and range between 10⁻² and 10⁻¹ M_⊙. 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 ∼5 × 10⁻² M_⊙ and ∼10⁻¹ M_⊙. 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 ∼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 ∼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 ∼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 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.

In this work, we presented the first collapse calculations of massive 1000 M_⊙ 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.

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_ϕ > 2v_r ,
• 2.  
  it should rotate faster than it falls vertically v_ϕ > 2v_z ,
• 3.  
  it should not be thermally supported $\tfrac{1}{2}\rho {v}_{\phi }^{2}\gt 2{P}_{\mathrm{th}}$, where P_th is the thermal pressure,
• 4.  
  it should be composed of dense material n > n_thre = 10⁹ cm⁻³, 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_system = 50 au from the analyzed sink. Note that we only analyze disks around primary stars (the most massive star within r < r_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_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₉₅ that encloses 95% of the mass. We also do the same in the vertical direction. Then we decompose the disk in N_θ = 50 regions around the polar axis in the mid-plane (within 4Δ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_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_thre = min(2 R_disk, R₉₅). The arbitrary pre-factor of 2 allows to preserve some of the non-axisymmetric features (e.g., spirals).

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

Zoom In Zoom Out Reset image size