© The Authors 2024

1 Introduction

Protostellar disks, often referred to as protoplanetary disks, are formed through the conservation of angular momentum during protostellar collapse. New observational evidences suggest that planets, or at least the gas giants, could form early during the evolution of those disks. The mass content of Class II–III disks indeed seems insufficient to explain observed exoplanetary systems (Manara et al. 2018; Tychoniec et al. 2020). In addition, the sub-structures of young < 1 Myr Class II (e.g., in HL-tau; ALMA Partnership 2015) and even < 0.5 Myr Class I (Segura-Cox et al. 2020) disks, particularly rings and gaps, could be indications of the presence of already formed giant planets. There are, of course, other theories for the formation of those structures (see the recent review by Bae et al. 2023), but the hypothesis of the presence of planets in gaps has recently been strengthened by kinematics evidence (Pinte et al. 2018, 2019). In contrast to older disks, Class 0-I disks could still have enough material to form planets. Unfortunately, the properties of these young disks are yet very poorly constrained. They are deeply embedded in a dense protostellar envelope, dominating the mass of protostellar objects during the whole Class 0 phase, and are often spatially unresolved in the wavelength range at which they can be observed (Maury et al. 2019; Sheehan et al. 2022).

From the theoretical perspective, we must resort to large scale simulations that self-consistently form disk populations. Significant efforts in the directon of such challenging modeling have been made by various teams in the pasts. For instance, Küffmeier et al. (2017), and later Küffmeier et al. (2019), investigated the impact of the accretion from large giant molecular cloud scales on the properties of disks, but without focusing on the formation of a full disk populations. This was initially done by Bate (2018), who investigated a full disk population forming in a massive protostellar clumps. Subsequently, Elsender & Bate (2021) investigated the impact of metallicity on disk population formation in very similar calculations, initially presented by Bate (2019). They mainly concluded that disk radii were decreasing with a decreasing metallicity. However, both studies did not account for the impact of the magnetic field.

Magnetic fields are, however, ubiquitous in observations of young stellar objects (YSOs; e.g., Girart et al. 2006; Rao et al. 2009; Maury et al. 2018). Observations suggest that they may play a key role in shaping the properties of some key features of the star formation process, such as the development of accretion flows, disk sizes and masses, and the occurrence of multiple stellar systems (Maury et al. 2018; Galametz et al. 2020; Cabedo et al. 2023). On the theoretical side, their role has been extensively investigated in the ideal (Price & Bate 2007; Mellon & Li 2008; Hennebelle & Fromang 2008; Hennebelle & Teyssier 2008; Joos et al. 2012) and non-ideal (Duffin & Pudritz 2009; Dapp & Basu 2010; Machida et al. 2011; Li et al. 2011, 2014; Dapp et al. 2012; Tomida et al. 2015; Tsukamoto et al. 2015; Marchand et al. 2016, 2020; Masson et al. 2016; Vaytet et al. 2018; Wurster & Bate 2019; Zhao et al. 2020, 2021; Hennebelle et al. 2020b; Mignon-Risse et al. 2021b,a) MHD frameworks for isolated collapse calculations of low- and high-mass cores. The magnetic fields have been proven critical to shape the disk through the regulation of angular momentum and for the launching of protostellar outflows. The importance of the large (clump-scale) magnetic field was also pointed out in the zoom-in simulations of Küffmeier et al. (2017, 2019). In the context of 50 M_⊙ mass clumps, Wurster et al. (2019), investigated the effect of all three non-ideal MHD effects on disk formation and concluded that they mostly impacted the small scales and the magnetic properties of the disks, but not their size and mass.

So far, only Lebreuilly et al. (2021) investigated disk formation in the MHD context with ambipolar diffusion for massive clump calculations, while systematically resolving the disk scales. These authors concluded that the clump scale magnetic field does indeed play a major role in setting the initial statistical conditions of the disks. Here, we build on this work by expanding the parameter space of the simulation suite with an overall higher numerical resolution. In this paper and the companion paper (Lebreuilly et al. 2024), we investigate in detail the impact of the initial clump conditions, namely, the magnetic field strength, treatment of the RT, and protostellar feedback (accretion luminosity, jets) and the clump mass. In this work, we will present six models and investigate the impact of the magnetic field and the RT modeling on the initial conditions of protostellar disks. This article is decomposed as follows. In Sect. 2, we briefly recall our methods, that are similar to those of Lebreuilly et al. (2021). In Sect. 3, we present in detail our fiducial model. This model will be our reference for comparison in this series of paper. In Sect. 4, we investigate the impact of the magnetic field and RT treatment on the initial conditions of our disk populations and their evolution. In Sect. 5, we describe the main caveats and prospects of our study, along with a first comparison with observation. Finally, we present our conclusions in Sect. 6.

2 Methods

2.1 Dynamical equations

To accurately describe the relevant physics in the context of star and disk formation, we solve the following dynamical equations: $$\begin{array}{c}\begin{array}{ll}\hfill & \frac{\partial \rho }{\partial t}+\nabla \cdot [\rho v]=0,\hfill \\ \hfill & \frac{\partial \rho v}{\partial t}+\nabla \cdot \left[\rho vv+\left(P_{\text{th}}+\frac{B^2}{2}\right)\mathbb{I}-BB\right]=-\rho \nabla \varphi -\lambda \nabla E_{\text{r}},\hfill \end{array}\\ \begin{array}{ll}\frac{\partial E}{\partial t}\hfill & +\nabla \cdot \left[v\left(P_{\text{th}}+E+\frac{B^2}{2}\right)-B(B\cdot v)\right]=-\rho v\cdot \nabla \varphi \hfill \\ \hfill & +{\text{Λ}}_{\text{AD}}-{\mathbb{P}}_{\text{r}}\nabla :v-\lambda v\nabla \cdot E_{\text{r}}\hfill \\ \hfill & \begin{array}{l}+\nabla \cdot \left(\frac{c\lambda }{\rho {\kappa }_{\text{R}}}\nabla E_{\text{r}}\right)+S_{\star },\\ \frac{\partial E_{\text{r}}}{\partial t}+\nabla \cdot \left[v{E}_{\text{r}}\right]=-{\mathbb{P}}_{\text{r}}\nabla :v+\nabla \cdot \left(\frac{c\lambda }{\rho {\kappa }_{\text{R}}}\nabla E_{\text{r}}\right)\\ \begin{array}{cc}& +{\kappa }_{\text{P}}\rho c\left(a_{\text{R}}{T}^4-E_{\text{r}}\right)+S_{\star },\\ \frac{\partial B}{\partial t}-\nabla \times [v\times B]& \\ & \begin{array}{l}+\nabla \times \frac{{\eta }_{\text{AD}}}{|B{|}^2}[((\nabla \times B)\times B)\times B]=0,\\ \begin{array}{cc}\nabla \cdot B& =0,\\ \text{Δ}\varphi & =4\pi \mathcal{G}\rho ,\end{array}\end{array}\end{array}\end{array}\hfill \end{array}\end{array}$$(1)

where ρ and v, E, E_r, and B are the gas density and velocity, total energy, radiative energy, and magnetic field. We also define the thermal pressure, P_th, gravitational potential, ϕ, radiative pressure, ℙ_r, Rosseland, κ_R, and Planck opacities, κ_P, as well as the radiative flux limiter, λ (Minerbo 1978), the temperature, T, the total luminosity source term, S_★, the ambipolar resistivity, η_AD, and Λ_AD is the heating term due to ambipolar diffusion. Finally, we define the gravitational constant 𝒢, the Stefan–Boltzmann constant a_R and the speed of light c.

To solve these equations, we use the adaptive mesh-refinement code (AMR) ramses (Teyssier 2002; Fromang et al. 2006) with RT in the flux-limited diffusion (FLD) approximation (Commerçon et al. 2011b, 2014), non-ideal MHD, and more particularly ambipolar diffusion (Masson et al. 2012) and sink particles (Bleuler & Teyssier 2014). More details on the code and modules that we used in this study can be found in the works referenced above.

2.2 Initial conditions

Our clumps are initially uniform spheres of 500–1000 M_⊙, with a temperature, T₀ = 10 K, and with an initial radius given by the thermal-to-gravitational energy ratio, α, such that: $$\alpha \equiv \frac{5}{2}\frac{R_0{k}_{\text{B}}{T}_0}{\mathcal{G}{M}_0{\mu }_{\text{g}}{m}_{\text{H}}},$$(2)

with k_B being the Boltzmann constant, m_H as the Hydrogen atom mass, and µ_g = 2.31 as the mean molecular weight. Owing to our choice for the values of α, all our clumps have the same initial radius of ~0.38 pc. The box size, L_box, is set to be four times larger, that is, L_box = 1.53 pc. Outside of the clump, the density is divided by 100. We set an initial turbulent velocity at Mach 7 with a Kolmogorov powerspectrum of k^−11/3 and random phases to mimic the molecular cloud turbulence. We point out, as explained in Lee & Hennebelle (2018), that the initial choice of spectrum for the turbulence has little impact on the results because the initial conditions are quickly forgotten has the collapse proceeds. Better ways to model the turbulence would require to start the simulation from even larger (kpc) scales, which is clearly beyond the scope of the present work.

The magnetic field strength is initialised according to the mass-to-flux over critical-mass-to-flux ratio µ such as: $$\mu =\left(\frac{M_0}{\varphi }\right)/{\left(\frac{M}{\varphi }\right)}_c,$$(3)

with ${\left(\frac{M}{\varphi }\right)}_c=\frac{0.53}{\pi }\sqrt{5/\mathcal{G}}$ (Mouschovias & Spitzer 1976).

2.3 Sink particles

Sink particles, following the implementation detailed in Bleuler & Teyssier (2014), are employed to mimic the behaviour of stars in our models. They are formed when the local density reaches the density threshold, n_thre = 10¹³ cm⁻³. This value is chosen in accordance with the analytical estimate of Hennebelle et al. (2020b). Once we form a sink, it automatically accretes, at each timestep, the material with a density above the threshold and within the sink accretion radius, 4∆x.

2.4 Radiative transfer modelling

We go on to examine two possible ways of modelling the RT. For all models, except one (NMHD-BARO-M500), we include the RT in the FLD approximation, using the solver of Commerçon et al. (2011b, 2014).

In the first approach, we considered sinks and stars as sources of luminosity, defined as the sum of the intrinsic and accretion luminosity. The former is computed using the evolutionary tracks of Kuiper & Yorke (2013), while the latter is defined as: $$L_{\text{acc}}=f_{\text{acc}}\frac{\mathcal{G}{M}_{\text{sink}}{\dot{M}}_{\text{sink}}}{R_{\star }},$$(4)

where R_★ is the star radius, also extracted from the tracks of Kuiper & Yorke (2013), M_sink and ${\dot{M}}_{\text{sink}}$ are its mass and mass accretion rate, and 0 < ƒ_acc < 1 is a dimensionless coefficient. Also, ƒ_acc corresponds to the amount of gravitational energy converted into radiation, in this work, we explore two values for ƒ_acc equal to 0.1 and 0.5. We refer to Lebreuilly et al. (2021) for an explanation on how the luminosity source terms are implemented in the code.

In the second approach (run NMHD-BARO-M500), we did not use the FLD approximation but instead assume a barotropic equation of state (EOS) to compute the temperature. This is, of course, an oversimplification, but it is interesting for two main reasons. First, barotropic EOS models are have lower temperatures than FLD calculations (Commerçon et al. 2010) and they allow for investigations of the effect of temperature on the disk formation and evolution. Second, because they are simpler than a full radiative tranfer modeling, these EOS are still widely used by the community.

For our barotropic EOS calculation, we assume (as in Marchand et al. 2016): $$T=T_0\sqrt{1+{\left(\frac{n}{n_1}\right)}^{0.8}}{\left[1+\left(\frac{n}{n_2}\right)\right]}^{-0.3},$$(5)

with n as the gas number density, n₁ = 10¹¹ cm⁻³, n₂ = 10¹⁶ cm⁻³, and T₀ = 10 K.

2.5 Refinement criterion

We use the AMR grid of ramses which allows us to locally refine the grid according to the local Jeans length. More specifically, we use a modified Jeans length such as: $${\tilde{\lambda }}_{\text{J}}=\{\begin{array}{ll}{\lambda }_{\text{J}}\hfill & \text{ if }n<{10}^9 {\text{cm}}^{-3},\hfill \\ \min \left({\lambda }_{\text{J}},{\lambda }_{\text{J}}\left(T_{\text{iso }}\right)\right)\hfill & \text{ otherwise. }\hfill \end{array}$$(6)

This modification is convenient for studying disk formation, especially in the presence of feedback since it is independent from the temperature at T > T_iso ≡ 300 K in the dense and heated regions. In all our models, we impose ten points per modified Jeans lengths within each cell to prevent artificial fragmentation (Truelove et al. 1997). The cell size is computed as a function of the refinement level ℓ as $$\Delta x=\frac{L_{\text{box}}}{2^{\ell }}.$$(7)

Our resolution always ranges from ~2460 au (~0.012 pc) in the coarsest cells of the simulation down to ~ 1.2 au (~5.8 × 10⁻⁶ pc) in the fine cells.

2.6 Disk selection

The disks analyzed in the present study were selected using the same method as in Lebreuilly et al. (2021), but we slightly modified the pre-selection criteria of Joos et al. (2012). As a reminder, the Joos criterion are: 1) n > 10⁹ cm⁻³, where n is the number density; 2) υ_ϕ > 2υ_r, υ_ϕ > 2υ_ɀ, where υ_r, υ_ɀ, and υ_ϕ are the radial, vertical, and azimuthal velocities, with the rotation axis being the direction of the angular momentum at the sink vicinity; and 3) $1/2\rho {\upsilon }_{\varphi }^2>2P_{\text{th}}$, where ρ is the gas density and P_th is the thermal pressure. In this work, we consider only the two first criteria. We have found that the last energy criterion arbitrarily removes the inner hot regions of the disks. Removing this criterion also allows a better comparison of models with a different accretion luminosity efficiency (and, hence, different temperatures).

Once all the disks of a model are selected, we analyze various of their internal properties (their radius, mass, temperature, etc.). For any quantity A, we compute the volume average such that: $$\langle A\rangle =\frac{{\displaystyle \sum _j{A}_j}\text{Δ}{x}_j^3}{{\displaystyle \sum _j\text{Δ}}{x}_j^3},$$(8)

where refers to all the disk cells of size, ∆x_j. The treatment of the temperature is slightly different, as we selected only the mid-plane cells to compute its averaged value. This allows for a better estimate of the temperature in the hot regions of the disk. In the remaining of the manuscript, we drop the 〈〉 notation foraverages as no confusion with the local value is possible. In addition, we estimate the disk radius as the median of the maximal extent in 50 equal-size azimuthal slices and the disk mass as the sum of the mass of every disk cell.

2.7 List of models

Our full list of models computed for this work is presented in Table 1. From left to right the table shows the model name, the initial clump mass, the thermal-to-gravitational energy ratio α, the mass-to-flux ratio µ (★ means ideal MHD), the accretion luminosity efficiency ƒ_acc (if applicable), the choice of RT modeling, the final median sink mass and the final SFE and corresponding time t_end.

3 Presentation of the fiducial model

Our reference model NMHD-F01 is a 1000 M_⊙ clump with µ = 10, ƒ_acc = 0.1 and a Mach number of 7. This run is essentially the same as the fiducial calculation of Lebreuilly et al. (2021), but with our new improved refinement criterion.

3.1 Large scales and star formation

We begin the description of our fiducial run NMHD-F01 by briefly presenting its evolution at the global scale. Figures 1a–c show the column density snapshots at various evolutionary stages for this model. The sinks are represented by the star markers.

As expected, the gravoturbulent motions lead to the formation of a network of highly non-homogeneous filament structures along which star formation mainly occurs (as in other similar studies, e.g., Lee & Hennebelle 2018; Bate 2018; Lebreuilly et al. 2021; Grudić et al. 2021; Lane et al. 2022). It is also very clear that the stars do not form in isolation here. Star formation is, in fact, more concentrated around one compact hub in the bottom-left part of the clump. This effect, which is most certainly a consequence of the global collapse, is off-centered because of the turbulence, which was also observed in the models of Lebreuilly et al. (2021) and Hennebelle et al. (2022). The global collapse of the clump is quite noticeable and non-isotropic, as expected in the presence of turbulence, which explains the presence of a main star-forming filament. This filament, connected with the previously mentioned hub, is the second-most active site of star formation of the clump.

We go on to focus on the main hub evolution. Figures 1d–f show the column density of NMHD-F01 at the same SFE as the top panels, but centered around the hub and at a smaller scale (12.5% of the box). Even at those scales, stars are clearly formed in filamentary structures which are connected to the larger scale network seen in the top panels of Fig. 1. These filaments are similar to the bridge structures that were observed in Küffmeier et al. (2019). They connect sinks with their neighbours and represent a shared reservoir of mass. They are relatively quiescent and typically survive a few ≃ 10 kyr. Quite clearly, we can see that a compact and highly interacting protostellar cluster is formed at the center of the hub. We point out that although sinks can get quite close to each other in this hub, we chose never to merge the sinks in our models since we are not resolving the scales for the stellar radii. This hub is a favored place for forming massive stars in the clump. In fact, the most massive stars formed in the model are part of this cluster.

Between SFE = 0.015 (t = 97.2 kyr) and SFE = 0.15 (t = 116.6 kyr), we observe a clear thickening of the filaments due to the radiative feedback of stars that heat up the gas and therefore increase its thermal support over time. This increase of thermal support significantly reduces fragmentation and sink formation, which essentially halts after a very efficient early phase (Hennebelle et al. 2022). Over the course of the simulation, integrated up to SFE = 0.15, about 90 sinks are formed, half of which are either single star or primaries (according to the simplified definition of Lebreuilly et al. 2021).

We point out that this model has formed fewer stars than its lower resolution counterpart (see the nmhd model of Lebreuilly et al. 2021). Very interestingly, the overall higher resolution of NMHD-F01 allows for the formation of one massive ~15 M_⊙ star over the course of the simulation. This is more massive by a factor of a few than the ones obtained in lower resolution runs (Hennebelle et al. 2022). We stress a clear correlation between sink masses and the mass of their surrounding envelope at a 1000 au scale, which indicates that the more massive star form in the more massive environment (see also, Klessen & Burkert 2000; Colman & Teyssier 2020). We note that for more extended dedicated descriptions (temperature, magnetic field and stellar mass spectrum) of the clump scales and star formation in very similar calculations, we refer to Hennebelle et al. (2022) and references therein.

3.2 Disks and small scales

In the following, we further describe the formation of structures, their properties, and their evolution at disk scales. There is a clear variety of disks and a wide diversity of commonly observed small scale features in the NMHD-F01 model and, more generally, in all of our models. Here, we describe the typical appearance and structures of our disks. As a support for that description, we show in Fig. 2 the edge-on or mid-plane density slices of 12 of our fiducial model disks at various times.

Compact disks. We observed many compact disks (see panels a–d, f, i, j, k, and 1). In fact, they represent the majority of the disks in the model, as we show later in this section (half of the disks are smaller than ~28 au at birth). As explained in Lebreuilly et al. (2021), they are a clear consequence of the regulation of the angular momentum by the magnetic braking during the protostellar collapse. It is worth mentioning that disks are indeed expected, from the observations, to be compact at the class 0 and early class I stage (Maury et al. 2019). The frequent occurrence of these < 50 au disks is compatible with the self-regulated scenario of Hennebelle et al. (2016). In this scenario, it is expected that the interplay between magnetic braking and ambipolar diffusion mostly leads to the formation of compact disks.

Substructures. Spirals and arcs (see panels g, h, and k) are also often observed, particularly in the presence of multiple systems and/or when the disk is gravitationally unstable and fragmenting. The latter effect is however restricted to the most massive and young disks (while they are still cold) as it is quite efficiently suppressed by the stellar feedback. Noticeably, ring structures are not present in our models: they are often attributed to planets-disk interactions and could, in principle, form from MHD instabilities in the disks (see a recent review by Lesur et al. 2022). The fact that we do not observe them may be attributed to two issues, namely: that either our resolution is still insufficient for them to occur in our disks or the early conditions in the disks (i.e., hot disks, with a massive turbulent envelope) are not favorable for rings to form. In general, the disk sub-structures are quite faint unless they originate from multiplicity in the disk (e.g., panel h). This is most likely a consequence of the thermal support due to radiation that stabilises the disk structure.

Magnetised flows. A common consequence of the magnetic field interplay with the gas at the disk scale is the triggering of interchange instabilities in some cases (as revealed by the prominent loop seen in panel c). This instability, which can transport momentum away from the disks (Krasnopolsky et al. 2012), was also observed in the rezoomed models of Küffmeier et al. (2017) in ideal MHD. Our study confirms that they are not always suppressed by the diffusive effect of ambipolar diffusion. Quite noticeably, and as also noted in Lebreuilly et al. (2021), magneto-centrifugal outflows and jets are absent in the models. This is most likely due to a suppression by the turbulence for the former (in accordance with Mignon-Risse et al. 2021a) and a lack of resolution in the inner regions of the disks for the latter. In Paper II, we will investigate the impact of protostellar jets on the disk properties using the sub-grid modeling of Verliat et al. (2022).

Non-axisymmetric envelopes. The protostellar envelopes around the disks are still massive by the end of the computation. Generally speaking, they are very diverse in structure and never spherically symmetrical, as previously pointed out by Küffmeier et al. (2017). From as far as a few ~1000 au scale up to the disk scales, accretion streamers are very common around disks in all the simulations. Those channels for high-density material accretion are typically (but not exclusively) connected to the disk mid-plane (e.g., panels a and b).

Flybys. Flybys, namely, close interactions between stars (not orbiting each other) or between a star and dense clumpy gas, are common during the clump evolution (Cuello et al. 2023, and references therein). Close flybys of two (or more) disks (panels i–k) are occurring several times in the model. In our simulations, they quite systematically lead to disk mergers, probably through the bridge structures (Küffmeier et al. 2019). Flybys have been shown to be able to truncate disks and trigger spiral formation in idealized calculations (Cuello et al. 2019). It is, however, difficult to establish their role in complex clump simulations that include many potential mechanisms to generate structures in the disks.

Disks in columns. Disks in columns, which are a very peculiar structure, with several occurrences at the early stages of the calculations (in all the models) only, is the formation of several disks in a single filament or column (panel 1). It is the consequence of the fragmentation of the clump on the filament scale. This is only possible when the latter are still cold. Indeed, this behaviour is not observed later when the stellar feedback heats up the gas and the thermal support precludes further fragmentation. The very short timescale and early disappearance (within less than 10 kyr after the first sink is formed) of these structures probably explains why they have not been observed.

Fig. 1 Evolution of the clump of run NMHD-F01 at various times (SFE = 0.0017, 0.015 and 0.15). Panels a–c: full column density maps in the (x–y) plane. Panels d–f: same maps but centered around sink 1 (located in the hub) and with an extent of 12.5% of the box scale. Sinks are represented by star symbols.

Fig. 2 Collection of disks from run NMHD-F01; edge-on or mid-plane density slices. In addition to the density, for each disk, we display the sink index, time of the corresponding snapshot, mass of the sink and of the disk, and the disk radius. Some disks, displayed at different times, may appear in several panels at once.

Fig. 3 CDF of the disks at birth time, as well as 10 and 20 kyr after for the NMHD-F01 model.

4 Disk populations

In the following, we discuss the properties of the disk populations for all of our models. We note that they were individually extracted for each model, as explained in Sect. 2.6.

4.1 Fiducial model

The disk populations of the NMHD-F01 model are described in detail in this section. In Fig. 3, we show the cumulative density function (CDF) of several key disk properties for run NMHD-F01 at the closest output from their birth time, but also 10 and 20 kyr afterward. We show the CDF of the radius (a), the mass (b), the ratio between the disk mass and the stellar mass (hereafter disk-to-stellar mass ratio, panel c), the mid-plane temperature (d), the magnetic field (e), and, finally, the plasma beta β ≡ 8/πP_th |B|² (i.e., the ratio between the thermal pressure and the magnetic pressure, panel f). Before describing those quantities individually, we point out that the disk populations are not strongly varying over time in the statistical sense except during the initial 10 kyr (hereafter the disk build-up phase). This does not mean however that disks are not individually evolving. As shown above, they are indeed non-linearly affected by interactions with the clump and other disks. In addition from the histograms, we show in Table 2 the mean, median (med.), and standard deviation (Stdev.) of these disk properties for NMHD-F01 at the birth time, at ages of 10 and also 20 kyr, and beyond.

It is worth mentioning that about ~25 disks are steadily detected in our fiducial model. This number, although it is lower than the findings of Lebreuilly et al. (2021), is consistent with the reduced fragmentation with our improved refinement criterion. The overall disk-to-star number ratio is higher than in this previous study as we now have ~70% of the systems hosting a disk.

4.1.1 Sizes and masses

Size. We first focus on the disk radii. This quantity is very interesting because it is perhaps the most reliable observable at all evolutionary stages. We can see in the radius CDF that the disks of NMHD-F01 are typically compact, with half of them being smaller than 30 au. This is in good agreement with the observations (Maury et al. 2019; Sheehan et al. 2022) and is slightly lower than what we found in Lebreuilly et al. (2021). This was expected, however, as we speculated that a higher resolution was needed to model the smaller disks. An interesting aspect of the evolution of the disk radius is that the radius of the smallest disk in the sample shifts from less than 10 au at birth to about 25 au for older disks. At the same time, the median value shifts from ~28 au to ~47 au. However, from 10 to 20 kyr, the PDF does not evolve significantly. To summarize, there is an initial build-up phase in the first 10 kyr of the disks lifetime during which they evolve a lot from compact to more extended disks as they (and their host star) accrete material from the envelope (Hennebelle et al. 2020b). Since this timescale is short, it would be almost impossible (for statistical reasons) to observe these disks. After this phase, the disk size does not change much over time, except in case of external perturbations (e.g., flybys or mergers) or when the system is a multiple.

Mass. It is of great importance to quantify accurately the disk mass since the mass content of the disks gives valuable information about the budget that is available for planet formation. An useful quantity to keep in mind is the minimum solar mass nebula (Hayashi 1981, MSMN). This mass, on the order of 0.01 M_⊙ and revised downward in more recent studies (e.g., Desch 2007), gives the minimal content needed to formsolar-like systems. In addition, the disk-to-stellar mass ratio provides an insight into the importance of self-gravity in the disk dynamics. Similarly to the disk radius distribution, there is a clear evolution of the disk mass and disk-to-stellar mass distribution during the build-up of the disks, namely, over ~ 10 kyr. During the early stages of their evolution, a significant fraction of mass in the system still belongs to the disk component which is comparable to (if not larger than) the mass of the stellar component. This is the stage during which the most massive disks can be gravitationally unstable, which could have consequences for early planet formation. After the bulk of the disk mass has been accreted by the protostar, namely, after a few kyr, the disk typically represents between 1 and 10% of the system mass. After this main build-up event, the disk masses are still typically larger than their initial value as the disk gets new material from the envelope. At this stage, half of the disks have masses between 0.01 M_⊙ and slightly less than 0.1 M_⊙, while the other half can reach masses up to ~0.3 M_⊙, which is still more than enough for planet formation according to the MSMN criterion. This value is, indeed, close to, if not below, the low disk mass limit of NMHD-F01 after the build-up phase, the mass of the disks in NMHD-F01 is most likely sufficient to form planetary systems similar to the solar system. A noteworthy detail is that the disk masses that we report are in good agreement with those of the hydrodynamical model of Bate (2018). Noticeably, the disk mass quite clearly correlates with the disk radius. This was of course expected as large volumes with the same typical density contain more mass. In Fig. 4, we show this correlation for the NMHD-F01 model. Each disk is displayed every 1 kyr and the various markers represent the different evolutionary stages of the disks (from birth up to 40 kyr for the older disks). The correlation between disk mass and disk radius is typically in between ∝ R_disk (plain line) and ∝ $R_{\text{disk}}^2$ (dotted line) and closer to the latter, which is expected for a perfectly symmetric disk. Figure 5 shows analogous information to Fig. 4 for the disk-to-stellar mass ratio versus the stellar mass. A clear correlation between the disk-to-stellar mass ratio and the sink mass is observed. It is slightly steeper, albeit close, to an inversely linear relation. We point out that this correlation also means that the disk mass is only weakly dependent on the stellar mass. In addition, as can be seen, the disks are more massive than the stellar component almost exclusively in the presence of low mass (below 0.1 M_⊙) and young (< l0kyr) systems. This is at odds with the previous study of Bate (2018), who has found a linear relationship between the disk mass and the stellar mass. We point out, as discussed extensively in Hennebelle et al. (2020b; see also Sect. 5.3), that the disk mass depends on the recipe used for the sink particles. This might explain why we did not observe the same correlation as found in Bate (2018).

Fig. 4 Disk radius versus disk mass for the NMHD-F01 model. Each disk is displayed once per kyr and the different markers and colors represent the various evolutionary stages for the disks.

Fig. 5 Same as Fig. 4, but for the disk-to-stellar mass ratio versus the primary mass. The horizontal dashed line represents a mass ratio of 1.

4.1.2 Temperatures

We now turn our attention to the mid-plane temperature distribution. The disks of NMHD-F01 are typically warm, with a median temperature of about ~100 K during the build up phase and about ~300 K at later stages. As shown in Hennebelle et al. (2020a) and later in Lee et al. (2021a), the temperature at the vicinity of a star can be controlled by its luminosity when it is sufficiently strong. In this context, we expect the disk temperature to scale as ∝ $L_{\text{acc}}^{1/4}$. The information that we show in Fig. 6 is analogous what is shown in Fig. 4, but for the correlation between disk temperature and their sink accretion luminosity. Those two quantities clearly have a correlation that is close to $T\propto L_{\text{acc}}^{1/4}$ (solid black line). We thus conclude that the accretion luminosity is indeed the dominant factor controlling the disk temperature in the model. More generally, we have confirmed this behavior for all the models, except NMHD-BARO-M500, for which the temperature is prescribed by the barotropic EOS.

It is worth mentioning that despite the weak correlation between the temperature and the accretion luminosity, it leads to a broad distribution of disk mid-plane temperatures (ranging from approximately 60 K to about 1000 K) since the accretion luminosity varies over four orders of magnitude. We have verified that the higher accretion rates and hence accretion luminosity correspond to the most massive stars of the system. As a consequence, the hotter disks are also those around the more massive stars in the model.

Fig. 6 same as Fig. 4 but for the mid plane temperature.

4.1.3 Magnetic fields

Finally, we bring our attention to the magnetic field strength and topology in the disks. First of all, half of the disks have a magnetic field below ~0.03 G during the build-up phase and below ~0.1 G after. We observe a slight increase of the magnetic field strength during the build-up phase which is probably due to an amplification through the still significant infalling motions that bend the field lines. The distribution of magnetic field is not varying significantly afterward. We also point out that none of the disks have a magnetic field greater than ~0.3 G. These are the two main consequences of ambipolar diffusion. The imperfect coupling between the neutrals and the ions indeed leads to a diffusion of the magnetic field at high density, where it is no longer dragged by the gas motions (see Masson et al. 2016; Mignon-Risse et al. 2021b,a; Commerçon et al. 2022, for similar observations in low- and high-mass isolated collapse calculations). A consequence of the low disk magnetisation is the negligible role played by the magnetic pressure inside the disk. The thermal pressure, enhanced by the stellar feedback from the accretion luminosity, is the main source of support (beside rotation) in the disk. We find that the magnetic pressure is typically about one order of magnitude lower than the thermal pressure and that a majority of disks exhibit 1 < β < 100. As we show later in this work, this is not the case for the ideal MHD model. It is noteworthy that (as for the magnetic field distribution) the distribution of β is not evolving much over time, except during the build-up phase.

In Fig. 7 (which is similar to Fig. 4 but for the disks vertical magnetic field vs. their radius), we clearly see an inverse correlation between the disk size and the vertical magnetic field strength, which very close to the ∝ |B_ɀ|^−4/9 (solid black line) predicted by the self-regulated scenario of Hennebelle et al. (2016). Despite being decoupled from the gas at the disk scales, the magnetic field does influence the disk size via the magnetic braking and perhaps via interchange instabilities outside of the disk, namely, in the envelope.

Fig. 7 same as Fig. 4 for the disk radius vs the absolute value of vertical magnetic field.

4.2 Impact of the magnetic field

We have computed two additional models with a different magnetic field treatment and strength. The first run, IMHD-F01, is the same as NMHD-F01 but evolving the magnetic field with an ideal MHD framework. The second run, NMHD-F01-MU50, is the same as NMHD-F01 but with a mass-to-flux ratio µ = 50, hence an initial magnetic field five times lower than the fiducial value.

In Fig. 8, we show the evolution of the disk quantities as a function of the SFE for NMHD-F01, IMHD-F01 and NMHD-F01-MU50. We show the disk radius (a), mass (b), mid-plane temperature (c), magnetic field strength (d), plasma beta (e), and the number of formed objects (stars, primaries, and disks, see panel f). For panels a to e, the dotted lines represent the median of the distribution and the colored surfaces represent the first and third quartiles of the distribution. Before describing each panel in detail, we note that the disk properties are still relatively steady with time and SFE for all three models (except during the build-up phase).

Fragmentation. We start with a focus on the number of objects, as shown in Fig. 8f. Quite clearly, the magnetic field actively stabilizes the cloud against fragmentation in NMHD-F01. Indeed NMHD-F01-MU50, which has an initially weaker magnetic field, is able to form more sinks (about 140, including about 60 primaries). Consequently, more disks are also formed in this run. Interestingly, almost all the systems are hosting a disk in the model, in contrast to NMHD-F01 , for which the fraction is closer to 75%. This is consistent with what was observed in Lebreuilly et al. (2021) where we have shown that the model without magnetic field (i.e., one with an infinite mass-to-flux ratio) forms more stars and more disks, and it has a higher fraction of disk-hosting stars. It is worth noting that star formation is happening more homogeneously in the NMHD-F01-MU50clump that is visibly more fragmented than that of NMHD-F01. This indicates that fragmentation is suppressed at a relatively large (larger than the disk) scale in NMHD-F01. This can be seen in Fig. 9, which shows the column density and the stars of NMHD-F01-MU50 at SFE=0.15; we can directly compare this with Fig. 1c. Similar observations were shown in Hennebelle et al. (2022).

If we now focus on the model with ideal MHD, IMHD-F01, we see that it forms even less stars than NMHD-F01. This is a consequence of the flux-freezing approximation that leads to an overall stronger pressure support of the magnetic field in the ideal MHD calculation. We notice that the fraction of primaries, which is the fraction of stars with no neighbours in a 50 au radius, and also the total number of primaries is, however, slightly higher that NMHD-F01 at the end of the calculation because small scale fragmentation is suppressed by stellar feedback.

Disk size. As we can see, while there are no strong differences of radii between the ideal MHD and the ambipolar diffusion case (as was also shown in Lebreuilly et al. 2021), a change in the mass-to-flux ratio however has an significant impact. The typical disk size is ~30–40 au in NMHD-F01 and IMHD-F01, and it is about twice as large in the case of a weaker initial magnetic field. In fact, disk sizes are actually slightly larger in IMHD-F01, which we explain later in this section. It is interesting to point out that this difference, by a factor of ~2, is perfectly consistent with the self-regulated scenario of Hennebelle et al. (2016); in their work, it is expected that the disk size is inversely scaled with the square root of the magnetic field strength. The present result is however a generalisation of what is proposed by Hennebelle et al. (2016) because, in the present case, the scaling also appears to be valid for the clump scale magnetic field. As pointed out by Seifried et al. (2013) as well as Wurster et al. (2019), the resemblance in terms of disk size between the ideal MHD and ambipolar diffusion models could be a due to turbulence. Indeed, it likely diminishes the influence of magnetic braking by generally producing a misalignment between the angular momentum and the magnetic field. However, as we have shown above, the clear correlation between the magnetic field strength and the disk size does indicate that they do play a crucial role in that regard. This is also supported by the fact that we obtain larger disk when considering a weaker magnetic field. We also emphasize that this similarity in the disk radius between the two models is actually misleading because the disk masses of IMHD-F01are lower than those of NMHD-F01.

Disk mass. The disks of NMHD-F01 -MU50 are more massive than those of NMHD-F01 and IMHD-F01 by a factor of about 2. This is not surprising since the disk masses are correlated to their radii in our models. Conversely, IMHD-F01 disks are generally less massive than those of NMHD-F01. The efficient magnetic braking at high density, in the absence of ambipolar diffusion, leads to more radial motions even in the disk, which allows for the star to accrete more material. This is a non-linear effect because more massive stars generate stronger feedback that also reduces fragmentation (Commerçon et al. 2011a). This aspect is interesting in the light of the similarity of the radii of the disk of NMHD-F01 and IMHD-F01. Essentially, even if their sizes are comparable, ideal MHD disks are typically less dense than the ones with ambipolar diffusion. The similarity in terms of disk size between ideal MHD and MHD with ambipolar diffusion is thus misleading and does not mean that magnetic braking is as efficient in the presence of ambipolar diffusion as it is without.

Disk temperature. Disks are hotter in IMHD-F01 than they are in NMHD-F01, the disks of NMHD-F01-mu50 being the coolest. The fragmentation in IMHD-F01 is suppressed, therefore, the stars are able to accrete more material which is further helped by a strong magnetic braking. Therefore IMHD-F01 stars have a stronger accretion luminosity, which leads to hotter disks. Conversely, fragmentation is most efficient for NMHD-F01-mu50 and, furthermore, the disks are larger, which means that a large part of their material is far away from the star and therefore colder.

Disk magnetic field and plasma beta. Here, we focus on the difference of magnetic field properties between the models. As a complement to Figs. 8d,e, here we show for the three models (from left to right): the correlation of the disk size versus the ratio between the vertical and azimuthal magnetic field (hereafter poloidal fraction, top panels) and the plasma beta versus the disk mass (bottom panels) in Fig. 10.

We see in Fig. 8d that at all SFE, the typical disk magnetic field is stronger in IMHD-F01 than it is in NMHD-F01 and unsurprisingly, it is the weakest in NMHD-F01-MU50. This explains quite clearly why the disks of NMHD-F01-MU50 are the largest in size. As we have shown previously, the disk size is roughly scaled as |B_z|−^4/9. We see in both NMHD-F01-MU50 and NMHD-F01 a weak correlation between the disk radius and the poloidal fraction, which is not clear at all for the IMHD-F01 model. This could be due to the more efficient winding-up of the magnetic field lines for the more massive disks, which is able to generate a significant toroidal field despite the diffusive effect of the ion-neutral friction that was observed in isolated collapse calculation of massive cores (Mignon-Risse et al. 2021b; Commerçon et al. 2022). As we can see in Fig. 10a, the bottom-right quadrant (small disks, strong poloidal fields) of the plot is dominated by the young (0–10 kyr) disks that did not have the time to wind up the field, while the top-left quadrant (large disks, weak poloidal fields) is dominated by older disks.

As mentioned earlier in this work, the previously described behaviour is not at all seen in the case of IMHD-F01; moreover, the poloidal fraction is generally lower in this model. This is a key difference that can be explained by the non-negligible impact of the ambipolar diffusion at the envelope scale for NMHD-F01. Without any diffusive effects, the magnetic field lines are already efficiently wound-up even before disk formation in the collapsing cores and, therefore, most disks already have a strong toroidal magnetic field at the birth time. Because of that, the correlation between the poloidal fraction and the disk size is not observed in the IMHD-F01 case. As explained earlier, the disk sizes of IMHD-F01 and NMHD-F01 are not strongly different. In fact (and perhaps counterintuitively), we find that ideal MHD disks are, on average, slightly larger than those of NMHD-F01. As explained in Lebreuilly et al. (2021), this is most likely a consequence of the strong toroidal fields, which stabilize the large disks against fragmentation and that are only formed in ideal MHD.

Quite interestingly,β is similar in the disks of NMHD-F01 and NMHD-F01-MU50, meaning that they reach the same level of magnetisation with respect to the thermal pressure. We clearly see from Fig. 10 that the beta plasma decreases with the sink mass for the three models; we also see that the typical value of beta is way lower for IMHD-F01 than it is for the two other models. The majority of IMHD-F01 disks are distinctly magnetically dominated (rather than thermally), whereas the opposite happens for NMHD-F01 and NMHD-F01-MU50. Low disk magnetisation is yet another clear consequence of ambipolar diffusion (see Masson et al. 2016; Mignon-Risse et al. 2021b,a; Commerçon et al. 2022, for similar results in low and high mass isolated collapse calculations). In ideal MHD disks, the infall actively drags the field lines, which causes a dramatic increase of the magnetic field intensity. In the case of IMHD-F01, because there is no diffusion at the envelope scale, the magnetic field is already strongly enhanced at the disk’s birth. The β-plasma then does not evolve much over time and stays, on average, close to 1, as can be seen in Fig. 8e. For the two other (MHD with ambipolar diffusion) models, the magnetic field is, as explained before, mostly vertical at birth. We then observed a slow decrease of β with time as a toroidal component is generated by the disk rotation.

We conclude that although the disk properties of the ideal MHD and MHD with ambipolar diffusion models present some similarities (i.e., radius, temperature), the difference in disk and stellar masses and magnetic field properties between the ideal MHD and MHD with ambipolar diffusion disks leads us to the conclusion that treating the ambipolar diffusion is crucial to better capturing the disk formation and evolution. We point out that knowing the initial configuration of the magnetic field is important for the onset and properties of MHD winds (see the review in Pascucci et al. 2023, and references therein). The strength of the magnetic field at the clump scale also appears to be an essential parameter that determines the properties of both the disks and star formation through its impact on the stellar feedback and magnetic pressure that acts against fragmentation.

Fig. 8 Evolution of the disk properties for NMHD-F01, IMHD-F01, and NMHD-F01-MU50 against the SFE. The lines correspond to the median of the distribution and the colored regions correspond to the area between the first and third quartile of the distribution. From left to right and top to bottom: radius, mass, mid-planet temperature, magnetic field strength, plasma beta, and number of objects.

Fig. 9 Same as Fig. 1c, but for the NMHD-F01-MU50 run. The clump is more fragmented as a result of the lower magnetic pressure support.

Fig. 10 Correlations between some disk and star properties for runs NMHD-F01, NMHD-F01-MU50, and IMHD-F01(from left to right). Top: disk size vs. ratio of poloidal over toroidal magnetic field. Bottom: beta plasma vs. primary mass. For the NMHD-F01-MU50and IMHD-F01panels, the information of the NMHD-F01 panel is duplicated with grey markers.

Fig. 11 Same as Fig. 8 for NMHD-F01, NMHD-F05, NMHD-F01-M500 , and NMHD-BARO-M500.

4.3 Impact of radiative transfer

In this section, we investigate the impact of the modelling of the temperature through the choice of the f_acc parameter and of the RT modeling (FLD, barotropic EOS). We ran the NMHD-F05 model, which is similar to NMHD-F01 but with f_acc = 0.5, therefore, we expected a more significant impact of the radiative feedback in this model. We also computed two additional models to better understand the impact of the RT modelling, NMHD-F01-M500 and NMHD-BARO-M500. To be comparable, both models have been run with α = 0.016, an initial clump mass of 500M_⊙, a Mach number of 7 and a mass to flux ratio of 10. However, contrary to NMHD-F01-M500, this has accounted for the RT with FLD and f_acc = 0.1, NMHD-BARO-M500 assumed a barotropic EOS. In a sense, NMHD-BARO-M500 gives an idea of what would be the evolution of the clump if no stellar feedback was included. We point out that this parameter choice also allows us to probe disk formation in less massive and less dense clumps with respect to NMHD-F01. In Fig. 11, we show the same information as in Fig. 8, but for NMHD-F01, NMHD-F05, NMHD-BARO-M500, and NMHD-F01-M500. In addition, the evolution of the accretion luminosity is shown in Fig. 12.

We first focus on the comparison between NMHD-F01 and NMHD-F05. If we look at panel f, we see that fragmentation is even more suppressed with f_acc = 0.5. We note that the effect of f_acc is non-linear: suppressing fragmentation leads to more massive stars that are brighter and hotter preventing fragmentation even more. The impact of f_acc is clear when looking at Fig. 8c, which shows the mid plane disk temperatures. The typical disk temperature quickly gets higher by a factor of ∼2 in NMHD-F05 than in NMHD-F01. The difference in temperature between the two models does not vary much over time which indicates that it is indeed caused by the accretion luminosity. Although f_acc is five times higher in NMHD-F05 than in NMHD-F01, the typical accretion luminosity of the former is almost one order of magnitude higher than in the latter. The additional factor of 2 in the accretion luminosity is mostly due to the reduced fragmentation in NMHD-F05 that leads to more massive stars.

Despite the difference among disk temperatures, other quantities such as the disk radius, mass, and magnetic field are similar in NMHD-F05 and NMHD-F01. This hints at the fact that the thermal support does not strongly affect the formation mechanism of the disks (magnetic braking vs. conservation of the angular momentum) and their evolution (quasi Keplerian rotation with a low viscosity as in an isolated collapse; Hennebelle et al. 2020b; Lee et al. 2021a). However, we note that the disks are typically thicker in NMHD-F05 than in NMHD-F01. This is not surprising at all since the disk scale height is controlled by the thermal support of the disk.

We now focus on the two 500 M_⊙ runs, NMHD-F01-M500 and NMHD-F01-M500. By the end of the calculation, at a SFE = 0.11, the barotropic EOS calculation has formed 140 sinks, whereas only 106 have been formed in NMHD-F01-M500. This is clearly an effect of the lower thermal support of the NMHD-BARO-M500, where the temperature is close to 10 K at low density. It is also interesting to point out that in the case of NMHD-BARO-M500, the maximum star mass is of ∼2 M_⊙, whereas it is around ∼10M_⊙ for NMHD-F01-M500. Despite being generally less massive (as disk fragmentation is more important in this model), the disks of NMHD-BARO-M500 are significantly more unstable due to their low thermal support. This can be seen in Fig. 13 that shows examples of fragmenting disks extracted from NMHD-BARO-M500. These disks are ubiquitous in the barotropic calculation, but rare for the other models with RT. This is an important point since disk fragmentation could be favourable for planet formation, either through streaming instability in the pressure bumps or through gravitational instability.

It is also worth mentioning that NMHD-F01-M500 and NMHD-BARO-M500 have a lower initial density that NMHD-F01. As a consequence, the sink accretion rates are generally lower and so is the accretion luminosity (see Sect. 5.2). If this has no consequences for NMHD-BARO-M500, which is computed with the barotropic EOS, the disks of NMHD-F01-M500 are impacted by this lower accretion rate and are colder than the ones of NMHD-F01. This confirms that controlling the accretion rate is crucial to predict the disk temperature and hence its fragmentation. As for the comparison between NMHD-F01 and NMHD-F05, it is interesting to point out that, temperature aside, the disk of NMHD-F01 and NMHD-F01-M500 are, however, quite similar. In addition for showing that the temperature is apparently not a controlling factor for the disk size, mass, and magnetisation (unless the disk are indeed very cold), this seems to indicate that the clump mass (or, rather, its initial density) is not a controlling factor either.

We conclude that a precise modeling of the RT including the impact of the accretion luminosity seems to be crucial to constrain the disks temperatures and the clump fragmentation. In addition, unless very cold disks are somehow relevant (if the barotropic EOS latter proves to be good approximation), the choice of accretion luminosity efficiency does not impact much the size and mass of the disks.

Fig. 12 Same as Fig. 11 for the accretion luminosity onto the disk-hosting primaries. The NMHD-F01, NMHD-F05 and NMHD-F01-M500 runs are displayed.

5 Discussion

5.1 Comparison with observations

One of the primary goals of the Synthetic Populations of Protoplanetary Disks project is to provide models for compare simulations and observations of Class 0 disks statistically.

We first present a tentative comparison of the disk populations with observed ones. For that comparison, we use a survey of disks in the Orion cloud (VANDAM survey; Sheehan et al. 2022). In Fig. 14, we show the cumulative distribution of the disk radius and mass of the populations extracted from all our models (coloured lines) compared with the observed ones (black lines). For the observations, we display the brute data of the survey with dashed lines and a version, re-scaled them by a factor 1/0.63, with plain lines. This re-scaling is done because the truncation radius as used in the (VANDAM survey, Sheehan et al. 2022) could be a better estimate of the radius containing 63% of the disk mass while our estimate better approximates the total disk radius. This crude re-scaling gives a good idea of what the total radius would be in this survey. In our populations, the disks are sampled every 1 kyr to mimic a diversity in evolutionary stages and to enhance the statistics in our clumps. We only selected the disks after the build-up phase, namely, after 10 kyr, to make sure that the distributions are in their steady phase, more likely to be observed. We also kept the populations of the different clumps separated to see the impact of the initial clump scale properties or physical assumptions on the disks.

In terms of disk radii, there is a moderately good agreement between our models and the observations. We note that the agreement is best for the models (with ambipolar diffusion) with a stronger initial μ = 10 magnetic field. Conversely, the NMHD-F01-MU50 model, which has an initial μ = 50, consistently produces discs larger than those observed. As was noted by Bate (2018), the radius that contains 63% of the mass could be a better measure of the radius in comparison with observations when a truncated power laws is assumed to fit the disks.

Considering that this value is more likely the one measured in the observations actually makes the agreement with the μ = 10 models even better. However, we point out that there is only a factor of ∼2 difference between the disk radius across all the populations. In addition it is important to stress a very important point. Disk observations are actually sensible to the continuum flux of the dust – and not the gas mass. The conversion from this flux to mass is not at all trivial, as we discuss below; therefore, it is not clear at all which of the r₆₃ or r_disk values (if any) are actually better probed by observations. This issue might be partly solved in Sheehan et al. (2022) though, as they that performed a careful radiative transfer modelling to fit the observed disks. We address this issue in a upcoming work, where we describe the post-processing of our models to produce synthetic observations. This will allow us to compare our models to real observations extracting the disks with the exact same methods.

We also recall that we found good agreement with disks radii from the CALYPSO survey (Maury et al. 2019) in the previous models of Lebreuilly et al. (2021). At that time, the magnetically regulated models were also in better agreement with the observations. This is also consistent with the observational evidence that show that only magnetically regulated models of the evolution of solar-type protostellar cores, with mass-to-fluxes ratio μ < 6 could reproduce the disk properties of the B335 protostar (Maury et al. 2018).

There is a more significant tension between our models and the observations when it comes to the disk mass. We generally report more massive disks than those obtained in Sheehan et al. (2022). We point out that similar tensions between models and observations were previously reported by Bate (2018, Fig. 25) and Tobin et al. (2020, Fig. 9) for the disk masses. As explained earlier, the masses obtained with our models are in line with those reported by Bate (2018) despite the significant differences in numerical methods. They indeed report, as we do, typical disk masses between < 0.01 and 1 M_⊙. We point out that the disks of Sheehan et al. (2022), from the Orion cloud actually have lower masses than those of Perseus (Tychoniec et al. 2020) and Taurus (Sheehan & Eisner 2017). Considering that the disk mass depends on the environment might partly solve the problem, as pointed out by Elsender & Bate (2021; their Fig. 9). At this stage, it is worth mentioning that the disk masses are poorly constrained both observationally and theoretically. Disk masses are likely controlled by the relatively unknown inner boundary condition namely, the star-disk interaction (see Sect. 5.3). On the observational side, the arbitrary choice of dust model and size distribution can lead to potentially large errors in the conversion between the flux from thermal dust emission measure at mm wavelengths, and the disks mass. This issue, discussed in a recent review by Tsukamoto et al. (2022), could be significant as the dust optical properties have been obsered to be very different in protostellar environments than in the diffuse ISM (e.g., Galametz et al. 2019). In addition, the computation of the gas mass from the dust relies on a conversion that assumes a constant gas-to-dust ratio, which is usually chosen to be 100. This hypothesis could be wrong in both ways: dust-rich disks (that can form if the dust decouples from the gas during the collapse Lebreuilly et al. 2020), where the gas-to-dust ratio could be lower than 100 and the disk mass would be overestimated, and dust-poor disks, for example, if some of the dust has already been converted into planetesimals, the gas-to-dust ratio would be larger than 100 and the disk mass would be underestimated. Unfortunately, molecular tracers might not lead to better estimates for the same reason as they also rely on a conversion factor to get the H2 disk mass. The gas kinematics could provide us with dynamical estimates of the disk mass assuming that the protostar mass is known (Reynolds et al. 2021; Lodato et al. 2023). Although this methods are challenging, they might be our best hope for a precise inference from the observational point of view.

It is also important to point out that most of the observed protostellar disks, including Class 0 disks, are older than 50 kyr, while our older disks are ‘only’ about 40 kyr by the end of our simulations. Their properties could, in principle evolve quite a lot through the class 0 objects, leading to an apparent disagreement between the mass estimates from models and observation that is, in fact, only due to an age difference. However, as we explained earlier, the disks properties are not significantly varying in our models for disks older than ∼10 kyr and, as was shown in isolated calculations (Hennebelle et al. 2020b), this relatively steady state probably lasts for at least 100 kyr. This supports us into thinking that there might indeed me a fundamental disagreement between models and observations for the disk masses.

Finally, and quite surprisingly, the NMHD-BARO-M500 model seems to be the model that actually fits the observations for both the mass and the radius much better. However, we stress that this likely is a coincidence. Barotropic models are indeed not supported by the recent e-disk survey (Ohashi et al. 2023) that have shown that young disks are quite hot and do not show obvious sub-structures except in rare and evolved case. This survey is thus is in much better agreement with our RT models.

Fig. 13 Collection of fragmenting disks from run NMHD-BARO-M500, with mid-plane density slices. In addition to the density, for each disk, we give the sink index, time of the corresponding snapshot, and mass of the sink and of the disk, as well as the disk radius.

Fig. 14 Comparison of the disk distributions of all our models (lines) with the observed ones in a sample of protostars in the Orion molecular cloud. A good agreement was found for the disk sizes but an important tension is observed when it comes to the disk mass. The observations uncertainties are not displayed for readability.

5.2 On the luminosity problem

The luminosity problem (Kenyon et al. 1990) is a long-standing issue of star formation. Observed YSO luminosities are below the values expected from steady-state protostellar accretion. This hints that either some of the accretion luminosity is not fully radiated away or the accretion is highly variable during protostar formation. Both of these issues can in principle be taken into account in the model through the efficiency factor, f_acc, which is a sub-grid modeling for the conversion of accretion luminosity into radiation in our calculations.

The typical accretion luminosities of our models with f_acc = 0.1 are typically a few 10 L_⊙, while they are about one order of magnitude higher in the case of NMHD-F05 that has f_acc = 0.5. Conversely, YSO observations seem to indicate lower luminosities with typical values on the order of a few L_⊙ during the Class 0 stage (Maury et al. 2011; Fischer et al. 2017).

As explained earlier, the accretion luminosity might not have a significant impact when it comes to the disk masses, sizes, and the properties of the magnetic field. However, constraining the accretion luminosity is of paramount importance for the disk temperature and our understanding of planet formation in those disks. First, the thermal support brought by stellar irradiation can act against the formation of structures in the disks (Rice et al. 2011, see also the comparison between NMHD-BARO-M500 and NMHD-F01-M500). In addition, the position of the snow lines depends on the temperature profile of the disk, which is important since planetesimal formation is expected to be most efficient in their vicinity, where the gas and dust properties of the disks abruptly change (see Drążkowska et al. 2023, and references therein). The position of the snowlines also determines the composition of the material available to form planets (gas and solids) at different locations in the disk. For instance, the location of the H₂O snowline is crucial to understanding in which conditions rocky planets form and how water is delivered to them. At the same time, in a relatively hot disk, some volatile species would never condense. Connected to that, a growing body of literature is studying how to link chemistry in disks to planet formation and to the composition of the exoplanets that we observe today (see for instance Öberg & Bergin 2021; Turrini et al. 2021; Pacetti et al. 2022). All of these works will considerably benefit from better constraints on f_acc.

At this stage, it is important to emphasize that the real and effective value of f_acc (if it is at all constant) could in principle be even lower than 0.1. In this case, there must either be an efficient mechanism to convert the gravitational energy into something else than radiative energy (e.g., magnetic energy or internal con-vective energy) or most of the radiation ought to be lost through the outflow cavity. As we show in Sect. 5.3, this issue could be tackled by models that resolve the star-disk connection up to small scales, namely, the stellar radii.

As previously mentioned, accretion in strong bursts could be invoked to solve the accretion luminosity problem (Offner & McKee 2011; Dunham & Vorobyov 2012; Meyer et al. 2022; Elbakyan et al. 2023). In fact, large variations of luminosity on timescales of years have been reported for some protostars (for example for B335; Evans et al. 2023). This is important because with a strong and steady accretion rate, the disk would be consistently kept warm; whereas short burst of accretion would not have long-term consequences on the disk temperature as the cooling timescale by the dust is very short (Hennebelle et al. 2020a) compared with the free-fall timescale. Figure 15 shows the mean accretion rate of the sinks as a function of time for NMHD-F01, NMHD-F05, NMHD-F01-MU50, and NMHD-F01-M500 (left) as well as the individual accretion rate of some sinks of NMHD-F01 (right). For all the models, there is a clear variability of the accretion rate over time. Accretion indeed occurs in short (but frequent) bursts of a few hundreds of years. During those bursts, the accretion rate rises to around 10⁻³M_⊙ yr⁻¹, occasionally reaching up to 10⁻² M_⊙ yr⁻¹. However, the average accretion rates are generally high in NMHD-F01 and NMHD-F05 in spite of these strong bursts; for these models, we observed typical average values around 10^−5–10⁻⁴M_⊙ yr⁻¹. This explains why, despite having a value of f_acc = 0.1, the accretion luminosity of our protostars is above typically observed values. Interestingly, the accretion rate is significantly lower in the case of NMHD-F01-M500 , where it is typically around 10⁻⁶–10⁻⁵ M_⊙ yr⁻¹. This clump being less dense than the fiducial run, it is also significantly colder (because it cools faster) and more efficiently fragmenting. Similarly, the weaker magnetic field of NMHD-F01-MU50 leads to more frag-mentation and lower accretion rates. We point out that episodic accretion is possibly (if not entirely likely) not fully resolved in our models both in terms of space and time. With shorter and stronger bursts, we could expect a lower typical average accretion rate, while still accreting enough stars to build the protostar.

It is not clear at this stage whether there is a luminosity problem in these models. The typical accretion rate indeed depends on the initial conditions of the clump and, as a result, the accretion luminosity can vary by a few orders of magnitude between the various models. We suspect that the various clumps that we present in this study are representative of different star-forming regions. With that in mind, it is worth pointing out that NMHD-F01/NMHD-F05 are probably more representative of massive star forming regions, whereas NMHD-F01-M500 is probably better at reproducing less compact nearby clumps. We do indeed have a top-heavy IMF in the case of NMHD-F01, contrary to the other runs, and as also shown in Hennebelle et al. (2022). An in-depth exploration of a larger parameter spaces in such simulations and of distant clumps in observations (e.g., Elia et al. 2017; Motte et al. 2022) would be required to confirm that hypothesis.

Fig. 15 Evolution of the stellar accretion rate as a function of time. Left: mean accretion rate for NMHD-F01 NMHD-F05, NMHD-F01-MU50 and NMHD-F01-M500 as a function of the time since the first star has formed. Right: same as the left panel, but for the individual accretion rates of a collection of sinks of NMHD-F01.

5.3 Caveats

5.3.1 The star–disk connection

A maximal resolution of ~1 au in the disks is already considerable with simulation boxes of ∼1.5 pc. Unfortunately, this minimum cell size still does not allow us to resolve the disk-star connection. To achieve that goal, we ought to be able to resolve the stellar radii, namely, to increase the resolution by at least two – if not three orders of magnitude. This is, of course, still impossible in the context of disk formation simulations. It is even more challenging in our case because we need to integrate the disks for a few tens of thousands of years to get a steady disk population. Because of that difficulty, we have to rely on the widely used sink particles (Bate et al. 1995; Federrath et al. 2010; Bleuler & Teyssier 2014) as they allow to integrate the model for a much longer time with a somewhat realistic inner boundary condition for the disk. It is important to recall that, unfortunately, the choice of the sink parameters (accretion threshold n_thre, sink radii and f_acc) does affect the calculation. Hennebelle et al. (2020b) showed that the mass of the disk was particularly affected by the choice of because the density at the center of the disk is quickly adjusting to this threshold. However, it is worth mentioning that the choice of this parameter is not completely arbitrary, we have indeed chosen n_thre based on the α-disk estimate of Hennebelle et al. (2020b). Fortunately, as they have shown, the disk radius is much less affected by the choice of n_thre probably because it is rather controlled by the magnetic field (Hennebelle et al. 2016).

In addition, we have shown in Sect. 4.3 that the impact of f_acc on the disk size, mass, and magnetic field is probably limited, provided that the accretion luminosity is not negligible (hence in all our models except for NMHD-BARO-M500). At the same time, f_acc does affect both the disk temperature and fragmentation.

To better constrain these two essential parameters, we strongly encourage studies dedicated to understanding the star-disk connection through the modelling of the stellar scales (Vaytet et al. 2018; Bhandare et al. 2020; Ahmad et al. 2023), while integrating the models in time as much as possible (Hennebelle et al. 2020b). This would provide the community with the necessary sub-grid modelling (star-disk connection) to constrain better the initial mass and temperature of protostellar disks.

5.3.2 Dust and planets

Recent studies indicate that dust evolution is not negligible during the protostellar collapse at disk-like densities (Guillet et al. 2020; Elbakyan et al. 2020; Tsukamoto et al. 2021; Bate 2022; Kawasaki et al. 2022; Lebreuilly et al. 2023; Marchand et al. 2023). This certainly carries important implications for the coupling between the gas and the magnetic field by means of the magnetic resistivities and probably on the RT because the opacities are dominated by the dust contribution.

In our collapse calculations, the dust size distribution is often assumed to be a constant in Mathis, Rumpl, Nordsieck (MRN) distribution (Mathis et al. 1977) and is only used to compute the resistivities and the opacity. Its evolution should, in principle, be taken into account self-consistently. This isstill very challenging in 3D simulations. In particular, in the context of large star-forming clumps, the memory cost of simulating multiple dust grains sizes would be extremely high. Fortunately, recent developments of new methods based on the assumption of growth purely by turbulence (Marchand et al. 2021, 2023) or accurate dust growth solvers that require only a few dust bins (Lombart & Laibe 2021) are paving the way to accounting for dust in future 3D simulations.

It is also worth mentioning that sufficiently large grains can, in principle, dynamically decouple from the gas (Bate & Lorén-Aguilar 2017; Lebreuilly et al. 2019, 2020; Koga et al. 2022). This phenomenon can lead to a local enhancement or depletion of the dust material, which would affect the dynamical back-reaction of the grains on the gas. In disks, this mechanism could also self-trigger the formation of substructures only in the dust (e.g., Dipierro et al. 2015; Cuello et al. 2019; Riols et al. 2020). Fully including the dust dynamics would also require to take into account the dust growth. This is clearly beyond the scope of this investigation.

Lastly, as our models do not follow dust evolution in the disks, we are thus unable to make predictions for planet formation and its implication. Again, this is beyond the scope of the present study. It is however important to keep in mind that fully formed planets or even planetary embryos could probably perturb the disk evolution and trigger the formation of substructures (e.g., Dipierro et al. 2016; Bae et al. 2017). To the best of our knowledge, planet population synthesis methods (see the review by Benz et al. 2014) are only used in 1D disks, however, their employment could be a way to tackle the problem in future works.

5.4 Hall effect and Ohmic dissipation

Ambipolar diffusion is not the only non-ideal MHD effect with a potential effect on disk formation processes. If the Ohmic dissipation is probably only dominant at very high densities (Marchand et al. 2016) and can thus be safely neglected, the Hall resistivity might be comparable, if not larger than the ambipolar resistivity in a significant density range in protostellar envelopes and possibly in the disks as well (Wurster 2021). For numerical reasons, we could not run our simulations with the Hall effect, but it is useful to recall its expected impact from our knowledge of isolated collapse calculations.

The Hall effect has indeed been investigated in detail by several groups over the past decade (e.g., Li et al. 2011; Tsukamoto et al. 2015, 2017; Wurster et al. 2016, 2019, 2021; Marchand et al. 2018, 2019; Wurster & Bate 2019; Zhao et al. 2020, 2021; Lee et al. 2021b). These works have typically found that the Hall effect could either enhance or decrease the rotation of the cloud (and the disk), depending of the initial relative orientation between the magnetic field and the angular momentum. Consequently, the Hall effect is expected to produce a bi-modality in disk properties. In the case where the Hall effect would enhance the rotation of the disk, therefore decreasing the effect of the magnetic braking, counter-rotating envelopes can typically be observed. It is worth noting that the Hall effect could also play a role in the fragmentation of the disks, when it accelerates rotation, as shown by Marchand et al. (2019).

It is not clear yet whether these effects are to be expected to play a strong role in the birth of disk populations obtained from a star-forming clump because of the dispersion in the relative orientation of the magnetic field and the angular momentum. So far, only Wurster et al. (2019) investigated disk forming in star-forming clumps with all three non-ideal MHD effects and have found no strong impact in the disk size. We point out that we also find a similar trend in our pure ambipolar diffusion case, although, as we pointed out, ideal MHD disk do have lower mass than the ones obtained with ambipolar diffusion. In addition, our ideal MHD disks do not have the same magnetic properties as the ones with ambipolar diffusion. In any case, simulations of massive star-forming clumps, including the Hall effect resolving the disks scales, would be very valuable for the community and should certainly be performed in the coming years.

6 Conclusion

In this work, we have explored the formation of protostellar disk populations in massive protostellar clumps with various assumptions and initial conditions. We now recall the main findings and conclusions of this work:

• Disk populations are ubiquitous in these simulations. A disk is found for around 70 to 90% of the stellar systems depending on the clump initial conditions;

• Disks are born across a variety of sizes, masses, structures, and nearby environments that reflect their individual history in a highly interacting gravo-turbulent star forming clump;

• We commonly find compact disks (in the presence of a strong magnetic field), non-axisymmetric envelopes (accretion streamers), substructures (spirals, arcs), magnetised flows (interchange instabilities), flybys, and peculiar structures such as disks formed in a single column. However, we find no ring structures or protostellar outflows or jets;

• Accretion luminosity is the dominant source of heating in the disks and it is controlling their temperature. The strength of the accretion luminosity depends on the clump properties. Clumps that fragment more efficiently also have lower accretion rates or luminosity, resulting in colder disks;

• The strength of the magnetic field at the clump scale is found to be a controlling factor for the disk size and the clump scale fragmentation. A stronger magnetic field leads to typically smaller disks and a reduced fragmentation. Because of this, the disks in clumps with a stronger initial magnetic field are also hotter;

• The accretion luminosity does not seem to be a controlling factor for the other disk properties (size, mass, and magnetic field). However, we point out that these properties could change if the disks were much colder than expected;

• The disk sizes obtained from our models are in relatively good agreement with the observed protostellar disk sizes from millimeter surveys. The case μ = 10 seems to be a better fit for the observations . However, there is tension with some surveys concerning the masses. Future post-processing of the models with radiative transfer tools should clarify the comparison between models and observations;

• Some well-known properties of isolated collapse calculations still hold in the context of large scale models. We confirm the important role of the magnetic field in shaping the disk masses and sizes and its combined importance with the radiative transfer in controlling their temperature and fragmentation. In addition, we show that when we account for ambipolar diffusion disks are weakly magnetized, while ideal MHD disks are not. Similarly to high mass star collapse calculations, we find that more massive disk generate stronger toroidal magnetic fields. Finally, we find that disks obtained in barotropic calculations tend to become fragmented more easily than those done using RT calculations. This confirms the validity of the isolated collapse approach to model protostellar disk formation, despite its incapacity to provide us statistics for the disk populations.

In this work, we have shown how diverse the populations of protoplanetary disks may be at their early stages and how they depend on their large-scale environments (magnetic field, radiation, cloud mass, etc.) as well as on the physical effects included (magnetic field with and without ambipolar diffusion, radiative transfer, etc.). We strongly encourage future works further exploring the influence of more clump properties (turbulence, size, and shape) as well as dedicated studies comparing such models with real observed data with synthetic observations produced with further elaborated radiative transfer codes.

Acknowledgements

We thank the referee for providing a very constructive report that helped us a lot to improve the quality of our paper. This research has received funding from the European Research Council synergy grant ECOGAL (Grant 855130). We acknowledge PRACE for awarding us access to the JUWELS supercomputer. This work was also granted access to HPC resources of CINES and CCRT under the allocation A0130407023 made by GENCI (Grand Equipement National de Calcul Intensif). G.R. acknowledges support from the Netherlands Organisation for Scientific Research (NWO, program number 016.Veni.192.233). Funded by the European Union (ERC, DiscEvol, project number 101039651). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. R.S.K. furthermore acknowledges financial support from the German Science Foundation (DFG) via the collaborative research center (SFB 881, Project-ID 138713538) “The Milky Way System” (subprojects A1, B1, B2 and B8), from the Heidelberg Cluster of Excellence “STRUCTURES” in the framework of Germany’s Excellence Strategy (grant EXC-2181/1, Project-ID 390900948), and from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 50OO2206). R.S.K. also thanks for computing resources provided by the Ministry of Science, Research and the Arts (MWK) of the State of Baden-Wurttemberg through bwHPC and DFG through grant INST 35/1134-1 FUGG and for data storage at SDS@hd through grant INST 35/1314-1 FUGG. M.G. acknowledges financial support from the French Agence Nationale de la Recherche (ANR) through the project COSMHIC (ANR-20-CE31-0009).

References

All Tables

All Figures

	Fig. 1 Evolution of the clump of run NMHD-F01 at various times (SFE = 0.0017, 0.015 and 0.15). Panels a–c: full column density maps in the (x–y) plane. Panels d–f: same maps but centered around sink 1 (located in the hub) and with an extent of 12.5% of the box scale. Sinks are represented by star symbols.
In the text

	Fig. 2 Collection of disks from run NMHD-F01; edge-on or mid-plane density slices. In addition to the density, for each disk, we display the sink index, time of the corresponding snapshot, mass of the sink and of the disk, and the disk radius. Some disks, displayed at different times, may appear in several panels at once.
In the text

	Fig. 3 CDF of the disks at birth time, as well as 10 and 20 kyr after for the NMHD-F01 model.
In the text

	Fig. 4 Disk radius versus disk mass for the NMHD-F01 model. Each disk is displayed once per kyr and the different markers and colors represent the various evolutionary stages for the disks.
In the text

	Fig. 5 Same as Fig. 4, but for the disk-to-stellar mass ratio versus the primary mass. The horizontal dashed line represents a mass ratio of 1.
In the text

	Fig. 6 same as Fig. 4 but for the mid plane temperature.
In the text

	Fig. 7 same as Fig. 4 for the disk radius vs the absolute value of vertical magnetic field.
In the text

	Fig. 8 Evolution of the disk properties for NMHD-F01, IMHD-F01, and NMHD-F01-MU50 against the SFE. The lines correspond to the median of the distribution and the colored regions correspond to the area between the first and third quartile of the distribution. From left to right and top to bottom: radius, mass, mid-planet temperature, magnetic field strength, plasma beta, and number of objects.
In the text

	Fig. 9 Same as Fig. 1c, but for the NMHD-F01-MU50 run. The clump is more fragmented as a result of the lower magnetic pressure support.
In the text

	Fig. 10 Correlations between some disk and star properties for runs NMHD-F01, NMHD-F01-MU50, and IMHD-F01(from left to right). Top: disk size vs. ratio of poloidal over toroidal magnetic field. Bottom: beta plasma vs. primary mass. For the NMHD-F01-MU50and IMHD-F01panels, the information of the NMHD-F01 panel is duplicated with grey markers.
In the text

	Fig. 11 Same as Fig. 8 for NMHD-F01, NMHD-F05, NMHD-F01-M500 , and NMHD-BARO-M500.
In the text

	Fig. 12 Same as Fig. 11 for the accretion luminosity onto the disk-hosting primaries. The NMHD-F01, NMHD-F05 and NMHD-F01-M500 runs are displayed.
In the text

	Fig. 13 Collection of fragmenting disks from run NMHD-BARO-M500, with mid-plane density slices. In addition to the density, for each disk, we give the sink index, time of the corresponding snapshot, and mass of the sink and of the disk, as well as the disk radius.
In the text

	Fig. 14 Comparison of the disk distributions of all our models (lines) with the observed ones in a sample of protostars in the Orion molecular cloud. A good agreement was found for the disk sizes but an important tension is observed when it comes to the disk mass. The observations uncertainties are not displayed for readability.
In the text

	Fig. 15 Evolution of the stellar accretion rate as a function of time. Left: mean accretion rate for NMHD-F01 NMHD-F05, NMHD-F01-MU50 and NMHD-F01-M500 as a function of the time since the first star has formed. Right: same as the left panel, but for the individual accretion rates of a collection of sinks of NMHD-F01.
In the text