© The Authors 2025

1. Introduction

In recent years, the role of magnetic fields in star formation has garnered a significant amount of interest. Aided by advances in far-infrared and submillimeter instruments capable of measuring linearly polarized dust emissions (e.g., ALMA, NOEMA, VLA), magnetic fields have been observed in dense cloud cores (Kirk et al. 2006; Jones et al. 2015; Kandori et al. 2018; Myers & Basu 2021) where they exhibit supercriticality (i.e., the mass-to-flux ratio is above unity) and a typical field strength of ∼10⁻⁵ G. Furthermore, using Zeeman line splitting techniques (Crutcher & Kemball 2019), they have also been observed in young stellar objects with values of ∼10³ G (Johns-Krull 2007; Johns-Krull et al. 2009; Yang & Johns-Krull 2011; Flores et al. 2024). Should the magnetic field be perfectly coupled to the fluid during the collapse of the dense core, which is the ideal magnetohydrodynamics (MHD) approximation, flux freezing implies that the resulting protostar would have a magnetic field strength of ∼10⁶ G, far in excess of observed values. Therefore, a considerable amount of magnetic flux is lost by the time the protostar becomes visible. This problem is known as the magnetic flux problem (Mouschovias 1985), which has so far eluded a concise answer. Current observational surveys of magnetic fields of young stellar objects (YSOs), although limited in sample size, have so far failed to find any correlation between magnetic field strength and stellar properties, such as their age and rotational period. However, they report a decreasing magnetic flux over time (Yang & Johns-Krull 2011).

The origin of the observed magnetic fields in YSOs is currently a subject of debate. Two main hypotheses dominate the discourse: the fossil field hypothesis, whereby the measured magnetic fields in these evolved sources are carried over from their inception in the second Larson core¹, and the dynamo hypothesis, which states that the measured fields are produced through a dynamo process. Ultimately, solving this problem requires a detailed model of the evolution of the protostellar magnetic field as the protostar transitions from the Class 0 to the Class I phase; the model should account for prestellar evolution and describe the magnetic field’s evolution using dynamo theory. However, such a model has not been developed yet, and little is reported on the subject in the literature. In the absence of any such model, the fossil field hypothesis remains favored and this may provide a constraint on star formation simulations as they must be able to form a protostar whose magnetic field has a strength of ∼10³ G.

A similar and closely linked issue is the angular momentum problem, which states that should angular momentum be conserved during the collapse of the dense core, stars would rotate far above their breakup velocity and circumstellar disks would be an order of magnitude larger than their observed sizes (∼30¹ AU, Maury et al. 2019; Tobin et al. 2020). Once again, the ideal MHD approximation fails to conform to observational data as it produces magnetic braking that is efficient enough to extract all angular momentum from the dense core (i.e., the magnetic braking catastrophe, Matsumoto & Tomisaka 2004; Hennebelle & Teyssier 2008; Hennebelle & Fromang 2008; Mellon & Li 2008). It is now widely admitted that resistive processes, mainly ambipolar diffusion, are responsible for breaking the ideal MHD limit toward higher density gas and reducing the magnetic braking efficiency to the point where a disk may form and reach sizes comparable to those found in observations (e.g., Masson et al. 2016; Hennebelle et al. 2016; Vaytet et al. 2018; Machida & Basu 2019; Wurster & Lewis 2020a,b; Mayer et al. 2024, also see the review by Tsukamoto et al. 2023). In addition, resistive MHD simulations report a converged magnetic field strength of ∼0.1 G in the first Larson core (primarily due to ambipolar diffusion), which allows the second Larson core to form with a magnetic field strength of ∼10³ G (e.g., Vaytet et al. 2018; Machida & Basu 2019; Wurster et al. 2022; Mayer et al. 2024).

In order to account for resistive processes, one must use a detailed chemical network that describes the abundance of charged species. In addition, one must also make an assumption on the dust grain size and density distribution in order to determine the surface area available for chemical reactions (Marchand et al. 2016, 2021, 2022; Zhao et al. 2020), and to account for the fact that the dust particles themselves may be the main charge carriers. In this regard, the Mathis-Rumple-Nordsiek distribution (MRN, Mathis et al. 1977), or some of its variants, is most often used as it is mostly valid for dust particles in the interstellar medium. However, recent studies having undertaken the effort of re-evaluating the dust size distribution during the collapse of dense cores have called into question the validity of the MRN distribution (Guillet et al. 2020; Silsbee et al. 2020; Lebreuilly et al. 2023; Kawasaki & Machida 2023; Tsukamoto et al. 2023; Vallucci-Goy et al. 2024). Most notably, these studies reveal an absence of small grains toward gas densities close to the first Larson core values (∼10⁻¹³ g cm⁻³), which in turn causes a stark drop in Ohmic resistivity, to the point where it is no longer a viable dissipative process at densities of the first Larson core and higher. The Hall effect, even within the MRN framework, remains the most poorly constrained resistive effect. Studies accounting for it report drastically different evolutionary scenarios for protoplanetary disks (e.g., Tsukamoto et al. 2015a; Wurster & Lewis 2020a; Wurster et al. 2022); however, the computed resistivities are too uncertain for us to draw any conclusions on the subject. The only resistive effect whose role and behavior during the protostellar collapse can be inferred with some confidence is ambipolar diffusion, whose reduction in magnetic braking efficiency allows the formation of disks whose sizes are in broad agreement with Class 0 disk size surveys (Hennebelle et al. 2016; Maury et al. 2019; Tobin et al. 2020). However, it still varies by orders of magnitude under different assumptions of dust coagulation.

Finally, it has recently become clear that subgrid models wrapped into sink particles that are used in protostellar collapse calculations in order to alleviate timestepping constraints produce results that are very sensitive to their parameters (Machida et al. 2014; Vorobyov et al. 2019; Hennebelle et al. 2020a). As these sink particles have a wide field of applications (e.g., Krumholz et al. 2009; Kuiper et al. 2010; Bate 2012, 2018; Hennebelle et al. 2020b; Mignon-Risse et al. 2021a, 2021b, 2023; Lebreuilly et al. 2021, 2024a,b; Commerçon et al. 2022; Grudić et al. 2022; Oliva & Kuiper 2023; Kuruwita et al. 2024), it is of vital interest to adequately constrain the subgrid parameters used in them. This requires one to study the innermost sub-au region of circumstellar disks, which entails second-collapse calculations resolving both the protostar and the newly formed circumstellar disk.

These second-collapse calculations can link the aforementioned issues together by providing a direct prediction of the magnetic flux and angular momentum inherited by the protostar, and by describing the star-disk interaction in detail. A number of these calculations exist in the literature (e.g., Banerjee & Pudritz 2006; Tomida et al. 2013; Vaytet et al. 2018; Machida & Basu 2019; Wurster & Lewis 2020b; Wurster et al. 2022; Mayer et al. 2024). However, most of these studies use idealized setups in which solid-body rotation is assumed and turbulence is absent in the initial dense cloud core. This absence of turbulence allows the magnetic field to maintain a coherence which amplifies its effects, be it magnetic braking or the launching of outflows and jets. In addition, with the exception of Machida & Basu (2019), Wurster & Lewis (2020b), Wurster et al. (2022) (who use nested-grid or smoothed particle hydrodynamics methods), these calculations are often stopped soon after protostellar birth owing to timestepping constraints, and the studies that have undertaken the challenge of integrating farther out in time have not provided an in-depth analysis of the gas kinematics in the innermost sub-AU region. Nevertheless, they found that the resulting protostar and circumstellar disk are thermally supported bodies, where thermal pressure gradient forces vastly outweigh their magnetic counterparts. In this regard, Ahmad et al. (2024) led a study in which the RHD approximation was used and found that the nascent protostar quickly reaches breakup speeds, at which point a circumstellar disk forms around it and expands outward. This occurs regardless of the initial conditions of the parent dense core. Machida & Matsumoto (2011) and Bhandare et al. (2024) similarly report the existence of this disk surrounding the protostar (hereafter called the inner disk or circumstellar disk).

Recent observational studies have also just begun probing deeply embedded Class 0 sources (Laos et al. 2021; Le Gouellec et al. 2024a,b). Although the structure of the system at the innermost sub-AU region is still difficult to infer from these observations, they seem to be reporting vigorous accretion onto the newly formed protostar.

In the present study, we expand upon our work in Ahmad et al. (2024) by including the effects of magnetic fields in our calculations, under the ideal and the non-ideal MHD approximations, while accounting for radiative transfer using the flux limited diffusion (FLD) approximation and including turbulence in our initial dense cloud core. Our goals are to describe the birth and early evolution of the protostar and the circumstellar disk surrounding it by focusing on the innermost sub-AU region. The simulations presented in this paper have the highest effective 3D resolution of all second-collapse RMHD simulations, all while pushing the calculations as far as possible in time. In light of the aforementioned recent studies around dust size distribution during protostellar collapses that report a stark drop in Ohmic resistivities, we chose to ignore Ohmic dissipation in our non-ideal MHD simulation, and only ambipolar diffusion is accounted for. We follow the collapse of a dense cloud core to stellar densities, and describe the initially isothermal phase of the collapse, the formation of the first Larson core and its subsequent adiabatic contraction, the second collapse following the dissociation of H₂ molecules, the birth of the protostar, and subsequently push the calculations as far as possible in time. In our pursuit of describing the smallest spatial scales relevant to protostellar and circumstellar disk birth, we obtained the best resolved protostars and circumstellar disk in the MHD literature. Particular attention is given to the structure of the magnetic field within the nascent protostar, and within the circumstellar disk. The evolution of the nascent circumstellar disk is also compared to its hydro counterpart in order to better ascertain the effects of magnetic fields on the system. Our results, reported below, carry multiple implications for the angular momentum and the magnetic flux problems. In addition, they offer constraints on subgrid parameters used in disk evolution studies.

In Sect. 2, we present our numerical and physical setup, whose results are discussed in Sect. 3. The implications of our results are discussed in Sect. 4. Finally, we conclude in Sect. 5

2. Model

2.1. Numerical setup

We used the RAMSES astrophysical code (Teyssier 2002) and its extension to MHD by Fromang et al. (2006). Non-ideal MHD (ambipolar + Ohmic dissipation) was implemented in the code by Masson et al. (2012), and radiative transfer under the flux limited diffusion approximation was implemented by Commerçon et al. (2011, 2014), González et al. (2015). The equation of state and opacity tables used were pieced together by Vaytet et al. (2013) from Saumon et al. (1995), Semenov et al. (2003), Ferguson et al. (2005), and Badnell et al. (2005). These account for a gas mixture consisting of 73% hydrogen and 27% helium, and assume a 1% dust mass content within the gas.

We use the same initial setup as in Ahmad et al. (2024), with the same refinement strategy. It consists of a uniform density dense cloud core of radius R₀ = 2.465 × 10³ AU, mass M₀ = 1 M_⊙ and temperature 10 K, equivalent to an initial ratio of thermal to gravitational energy of 0.25. Angular momentum is present through the inclusion of a turbulent velocity vector field parameterized by the turbulent Mach number ℳ, which we have set to 0.4. Radiative transfer is accounted for under the gray FLD approximation. A uniform magnetic field threads the dense cloud core along the z axis, and its strength is parameterized by the mass-to-flux ratio which we have set to 4. This corresponds to an initial magnetic field strength of ∼10⁻⁵ G in the dense cloud core, and an Alfvénic Mach number of ℳ_a ≈ 0.12. This setup is identical to that of run G2 in Ahmad et al. (2024), with the only difference being the presence of magnetic fields.

Two simulations will be presented in this section; one under the ideal MHD approximation (hereafter IMHD) and one in which we have accounted for ambipolar diffusion (hereafter NIMHD). Run IMHD is mainly presented in this paper for comparative purposes with the more realistic run NIMHD, as that allows us to better isolate the role of magnetic resistivities during the simulation.

These two simulations use the same refinement strategy; however, run IMHD has a maximum refinement level of ℓ_max = 26 (the coarsest level is at ℓ_min = 6) whereas run NIMHD has ℓ_max = 25. As we show later-on, this is because run IMHD forms a much more compact protostar, whose properties require a finer spatial resolution to describe. This means that at the finest refinement level, run IMHD and NIMHD respectively have a spatial resolution of Δx_IMHD = 1.4 × 10⁻⁴ AU and Δx_NIMHD = 2.9 × 10⁻⁴ AU.

2.2. Zoom-out

Owing to very stringent time-stepping constraints following protostellar birth, run NIMHD requires approximately two days of CPU wall time in order to integrate ≈40 hours. This is because the timestep at the finest level reaches a mere minute, and the poor load balancing causes most cells to be handled by a few CPUs. As the protostar and circumstellar disk grow and expand over time, this problem is aggravated as a considerable number of cells are created to describe the newly formed structures. In order to alleviate the timestepping constraints, we have also run a simulation branched out of run NIMHD nearly 0.4 years after protostellar birth, in which the maximum refinement level was reduced from ℓ_max = 25 to ℓ_max = 24. This allowed us to push the simulation considerably farther out in time², which is particularly useful to study the expansion of the newly-formed circumstellar disk. This run, labeled “NIMHD_LR”, is discussed in Sect. 3.4.

Run IMHD required 2 months and 4 days of CPU wall time, whereas run NIMHD required 7 months and 6 days. Run NIMHD_LR ran for 6 months and 13 days from its branch-out point, meaning that when including the zoom-out, the non-ideal simulation ran for over a year of CPU wall time. All simulations were run on 64 CPU cores.

3. Results

3.1. Large scale structures

We first begin by describing the system at the scale of the dense core itself (∼10³ AU), with the goal of providing the contextual environment in which the protostar is born. To this end, we compare runs IMHD and NIMHD in Fig. 1 at our final simulation snapshots (respectively ≈23.25 and ≈23.39 kyr after simulation start), which displays the column density (first row), the optical depth τ computed along the line of sight (second row), and the maximum temperature along the line of sight (third row). τ is computed as

$$\begin{array}{c}\hfill \tau ={\displaystyle {\int }_{z_{\min }}^{z_{\max }}}\rho {\kappa }_{\mathrm{R}}dz,\end{array}$$(1)

Fig. 1. Comparison of runs IMHD (first row) and NIMHD (second row) at the scale of the dense cloud core itself. The snapshots are taken respectively at t ≈ 23.25 and t ≈ 23.39 kyr following the collapse of the dense core. The first column displays column density (panels a and b), the second column displays the optical depth computed along the line of sight (panels c and d), and the last column displays the maximum temperature along the line of sight (panels e and f). All maps are projections along the z-axis. The green contour in panels (c) and (d) represent an optical depth of unity. The scale bars in the first column apply to the other columns as well.

where ρ is the gas density and κ_R the Rosseland mean opacity.

The column density maps show a filamentary structure of size ∼10² AU forming in both runs. This structure is formed by gravo-turbulence (Tsukamoto & Machida 2013); however, it appears much thinner in these calculations than their RHD counterparts in Ahmad et al. (2024). This is due to magnetic braking, which extracts a significant amount of angular momentum from the gas and thus prevents it from spreading out as much. Since ambipolar diffusion begins acting at higher densities (∼10⁻¹⁴ g cm⁻³), the two runs yield identical column density maps outside the filament; however, in the case of run NIMHD it has fragmented into two distinct dense cores (Fiege & Pudritz 2000). The existence of a secondary bound fragment within the filament is owing to an extended first core lifetime. The first core survived ≈100 years longer in run NIMHD owing to a reduced mass accretion rate onto it, which is in turn due to less efficient magnetic braking. In this time span, the filament fragmented in run NIMHD, whereas the stringent timestepping following the second collapse froze the simulation at larger scales in run IMHD, and no bound fragment is witnessed at its final simulation snapshot.

Despite the very similar structure, the two runs have differing optical depth maps (Figs. 1c and d). Run NIMHD has a more spatially extended optically thick region (lime-colored contours) than run IMHD. Since the two runs display identical column density values at the location where this optical thickness is achieved in run NIMHD, the differing optical depth maps are due to the differing temperatures found within the cloud core, as ambipolar diffusion significantly heats-up the gas (Fig. 1f). The increase in temperature at these densities manifests itself as an increase in opacity (see Figure 1 of Ahmad et al. 2023). This serves to show that the two models should produce distinct emission maps that may be discriminated against with current observational instruments.

We now turn to studying the collapse in quantitative terms using Fig. 2. Figure 2a displays the maximum density of the simulation as a function of time since first core formation (defined as the moment where ρ_max > 10⁻¹⁰ g cm⁻³). The steep rise in ρ_max in this figure corresponds to the second collapse (i.e., protostellar birth). We see here that the two runs display different first core lifetimes, with run NIMHD entering the second collapse phase nearly 200 years later. This is in contrast to Vaytet et al. (2018)’s results, who reported a longer first core lifetime in their ideal MHD simulation due to the interchange instability reducing mass accretion rates onto it. This discrepancy between our results is likely due to our use of turbulent initial conditions, which although does not prevent the emergence of the interchange instability in run IMHD (Fig. 1a), still reduces its efficiency³. The first core in run NIMHD survived for a total of ≈250 years, which is about half as much as the hydrodynamical run presented in Ahmad et al. (2024). Its extended lifetime in comparison to run IMHD is due to the reduced magnetic braking efficiency, which allows for angular momentum to reduce mass accretion rates onto the first core. The maximum density reached post-second collapse in run IMHD is ∼10⁻¹ g cm⁻³, and ∼10⁻³ g cm⁻³ in run NIMHD.

Fig. 2. Quantitative comparison of the collapse between run IMHD (black) and run NIMHD (red). Panel (a) displays the evolution of the maximum density since the formation of the first Larson core, which we define as the time when a density of ≈10−10 g cm−3 is achieved. Panels (b) and (c) display the magnetic field strength evolution as a function of time since first core formation and since protostellar birth (defined as the moment a density of ≈10−5 g cm−3 is reached), where the solid lines represent the maximum magnetic field strength and the dotted lines represent the field’s strength measured at the location of maximum density.

During the collapse, flux freezing causes the magnetic field strength to increase with increasing density (with B ∝ ρ^2/3). This causes the maximum magnetic field strength (B_max) in run IMHD, shown in Fig. 2b, to continuously increase over time (with increasing central density), with small drops in magnetic field strength being caused by turbulent reconnection. In run NIMHD, however, the magnetic field strength displays a plateau at ≈0.5 G, owing to ambipolar diffusion. This field strength within the first core has been consistently retrieved by numerous studies in the literature that account for ambipolar diffusion, in the low-mass and in the high-mass regime (e.g., Masson et al. 2016; Vaytet et al. 2018; Mignon-Risse et al. 2021b; Wurster et al. 2022; Mayer et al. 2024). Once the second collapse occurs, flux freezing (which is recovered in run NIMHD following dust sublimation and the ionization of atomic gas species) once again causes a strong increase in magnetic field strength, which reaches ∼10⁵ G in run IMHD and ∼10³ G in run NIMHD. We also notice in both runs that the magnetic field strength measured at the location of maximum density (B_central, dotted lines) is a factor of ≈2 below B_max. Following the second gravitational collapse (Fig. 2c), the maximum magnetic field strength reaches ≈3 × 10⁵ G in run IMHD, and ≈10⁴ G in run NIMHD. Soon after protostellar birth, the maximum field strength in run IMHD continuously decreases to ∼5 × 10⁴ G, and in the case of run NIMHD, it decreases to ∼6 × 10³ G and plateaus around this value. We see the same trend in B_central, which fails to coincide with B_max following the second collapse, and whose discrepancy with it seems to worsen over time. We show later in Sect. 3.3 that the drop in magnetic field strength in run IMHD is mostly due to an outward advection of magnetic flux. The discrepancy between B_max and B_central has been reported in previous papers in the literature, most notably Wurster & Lewis (2020b), Wurster et al. (2022). They also report a reduction in B_max shortly following protostellar birth. Our results confirm their findings; however, the cells containing ρ_max and B_max are separated by a very small distance (∼10⁻² AU) in our simulations, whereas in theirs it is of the order of ∼1 AU, and we do not report the existence of a “magnetic wall” on which magnetic flux is accumulated as they do.

3.2. The second collapse

3.2.1. Qualitative result of the second collapse

We now turn to the main focus of our study; the structure of the system following the second gravitational collapse. To this end, we first begin by studying the qualitative structure of the system with the aid of density, temperature, radiative flux, and radial velocity slices displayed in Figures 3 and 4. The slices are projected along the angular momentum vector of the gas within 0.2 AU for run NIMHD, and in the case of run IMHD, along the z axis since there is no angular momentum left in the gas owing to the efficiency of magnetic braking.

Fig. 3. Set of slices showing the evolution of the density (first two columns) and temperature (last two columns) for run IMHD (first and third columns, panels a–e and k–o) and run NIMHD (second and fourth columns, panels f–j and p–t). Each row represents a different time, where t = 0 corresponds to the moment of protostellar birth. For comparative purposes, the slices are shown at similar times, and the timestamp is written in the top right corner of each panel. The slices are done in the z direction for run IMHD, and along the angular momentum vector for run NIMHD. The scale bars in the first row apply to all the other rows as well.

Fig. 4. Same as Fig. 3, but now showing the radiative flux (first two columns) and radial velocity (last two columns).

The differences between the resulting protostars are stark. The first row displays the system at protostellar birth, which we define as our t = 0⁴. The protostar in run IMHD is more compact than run NIMHD, displaying higher densities and temperatures owing to the lack of centrifugal support against gravity. This causes it to form at a radius of ≈0.97 R_⊙, whereas in run NIMHD, the centrifugal support flattens the protostar considerably and extends its radius to ≈4.8 R_⊙. In the weeks following the formation of the protostar in run IMHD (panels b–e and l–o of Fig. 3 and Fig. 4), its size grows considerably as it accretes material from its surroundings. This is due to the subcritical nature of its accretion shock, which struggles to radiate the incoming accretion energy (see Ahmad et al. 2023). In addition, filamentary structures protruding from the stellar surface can be seen growing in spatial extent as time progresses. These are in fact current sheets akin to coronal mass loops, which appear as filaments when visualized in 2D slices.

In the case of run NIMHD, an entirely different evolutionary sequence is witnessed. As time progresses, a disk-like structure surrounding the protostar is formed. This is due to the latter’s accumulation of angular momentum, after-which it reaches breakup velocity and material is advected outward⁵. This result confirms the findings of Ahmad et al. (2024), in which no magnetic fields were present. As the disk grows in size and in mass, it exhibits spiral waves which form as a result of gravitational instabilities. These spiral waves carry a significant amount of angular momentum outward and cause increased mass accretion rates onto the protostar.

These figures show the importance of including ambipolar diffusion in our calculations, as it allows for enough angular momentum to survive and hence averts the magnetic braking catastrophe. Below, we proceed to providing a more quantitative analysis of the nascent structures, and since the two runs yielded drastically different results, we define the protostar in two different ways which are presented in Appendix A.

3.2.2. Gas structure and kinematics

Having ascertained the qualitative structure of the system in the two runs, we proceed to providing a more quantitative comparison between run IMHD and run NIMHD. To this end, we display in Fig. 5 and Fig. 6 averages of various physical quantities. In the case of run NIMHD, since the structure we witness is a flattened disk-like structure, these quantities are averaged azimuthaly in cylindrical bins in which only cells in the midplane region are selected. The midplane is defined as the region in which z ∈ [ − 2.5;  2.5]×10⁻² AU, where the z component is computed along the angular momentum axis of the gas within 0.1 AU. In the case of run IMHD, the measurements are done using the spherical coordinate system since the protostar and the distribution of material around it posses a spherical morphology.

Fig. 5. Set of measurements of various physical properties of run IMHD at protostellar birth (t = 0, solid lines) and t ≈ 117 days later (dotted lines). These are averages in spherical bins, and show the gas density (panel a), temperature (panel b), entropy (panel c), radial and azimuthal velocity (panels d and e), magnetic field intensity (panel g), and the azimuthal and meridional components of the magnetic field, normalized by its magnitude (panels h and i). Panel (f) displays the enclosed magnetic energy (black lines) and the enclosed mass (red lines) as a function of spherical radius.

Fig. 6. Set of measurements of various physical properties of run NIMHD. These were done in cylindrical radial bins, in which only cells belonging to the midplane, defined as z ∈ [ − 2.5;  2.5]×10−2 AU, were used. The solid lines are measurements done at t = 0 (corresponding to the moment of protostellar birth), and the dotted lines are measurements done ≈190 days later. Panels (a), (b), and (c) display respectively the gas density, temperature, and specific entropy. Panels (d) and (e) display the gas (cylindrical) radial and azimuthal velocity. The red curves in panel (e) display the Keplerian velocity computed with the protostar’s mass at a given snapshot. Panels (g), (h), and (i) display the magnetic field strength, its toroidal component, and its vertical component. The toroidal and vertical components are normalized by the total magnetic field strength. Panel (f) displays the enclosed magnetic energy (black lines) and the enclosed mass (red lines) as a function of spherical radius.

We begin by studying the structure of run IMHD (Fig. 5). The density profile at protostellar birth, displayed in Fig. 5a (solid line), shows that the central region of the protostar reaches ∼10⁻¹ g cm⁻³⁶. Nearly 117 days later (dotted line), this value drops to ∼10⁻³ g cm⁻³. In both snapshots, a power-law tail follows the central density peak. The sharp discontinuity in the radial velocity profile (Fig. 5d) displays the location of the accretion shock, which is ≈4 × 10⁻³ AU (≈0.86 R_⊙) and subsequently moves outward as the protostar expands. The azimuthal velocity curves, shown in Fig. 5e, display the efficiency of magnetic braking in this simulation: nearly no angular momentum survived, as v_ϕ alternates between positive and negative values and is mostly a noisy measurement. In Fig. 5c, the specific entropy of the gas⁷ is shown. Here, as in Ahmad et al. (2023), we see that ds/dr > 0 throughout the protostellar interior, meaning that the protostar is radiatively stable against convective instabilities. However, as accretion still drives turbulence within the protostellar interior (Bhandare et al. 2020; Ahmad et al. 2023), the entropy profile at our final snapshot is flattened as a result of the mechanical transport of energy. The entropy profile has also been lifted upward as a result of the accretion of energy. These results are very similar to the spherically symmetrical RHD run presented in Ahmad et al. (2023), which again illustrates how efficient magnetic braking is in this run.

In Fig. 6, however, we again witness very different results for run NIMHD. First, the density reached in the central regions is two orders of magnitude lower and at ∼10⁻³ g cm⁻³, a value close to the hydro runs presented in Ahmad et al. (2024). Unlike in run IMHD, however, no hydrostatic bounce occurs and the maximum density remains constant as time progresses. The temperature shown in Fig. 6b is also an order of magnitude lower than in run IMHD, and sits close to 5 × 10³ K; however, this increases to 7 × 10³ K nearly 191 days later, meaning that the protostar is heating up as it accretes material. The cylindrical radial velocity v_cyl, displayed in Fig. 6d, shows that the protostellar accretion shock is formed at 2 × 10⁻² AU (≈4.3 R_⊙), which is nearly five times larger than in run IMHD. This midplane shock front expands outward to 2 × 10⁻¹ AU at our final snapshot. We emphasize that it is no longer the protostellar accretion shock that is displayed by the discontinuity in v_cyl, but rather, that of the newly-formed circumstellar disk in which the protostar is embedded. In Fig. 6e, we display the azimuthal velocity curves. Here, we see very clearly that rotational motion exists, as v_ϕ > 0 throughout the radii displayed in the figure. Furthermore, these curves show that the central regions are in solid body rotation, whereas the circumstellar disk exhibits differential rotation. At t ≈ 190 days, the disk displays a fully Keplerian ($v_{\mathrm{K}}=\sqrt{G{M}_{\ast }/r}$) rotation profile (dashed red curve).

Finally, we display in Fig. 6c the specific entropy of the gas. As in run IMHD, the protostar is radiatively stable against convective motion. However, at our final simulation snapshot, we see that it is no longer the case as there exists a region in which ds/dr < 0. The existence of this region is likely due to the prominent spiral waves within the disk, and the negative entropy gradient means that a convective instability occurs (Schwarzschild 1906), which further contributes to turbulent motion within the disk and hence enhances accretion onto the protostar.

This analysis once again shows the stark differences of both runs; whereas the almost entirely complete absence of angular momentum in run IMHD owing to magnetic braking causes the second collapse to form structures more akin to those produced in spherically symmetrical calculations, the inclusion of ambipolar diffusion allows a considerable amount of angular momentum to survive and hence form a rotationally supported disk surrounding the protostar. In this sense, run NIMHD is more related to hydrodynamical runs than to run IMHD, and thus should be quantitatively and qualitatively compared as such.

3.3. Magnetic field structure

In this section we describe the structure and morphology of the magnetic field within and in the close vicinity of the protostar. To this end, we display in Fig. 7 and Fig. 8 slices showing the magnetic field strength and plasma β (= 8πP/B²). In Fig. 7, the slices are shown in a top-down (panels a–e and k–o) and edge-on (panels f–j and p–t) view, whereas the absence of rotation in run IMHD renders any such double-visualization useless, and we may visualize these quantities along the z-axis only. We also leverage the information available in Fig. 5 and Fig. 6.

Fig. 7. Magnetic field of run NIMHD: a set of slices displaying the magnetic field strength (top two rows) and plasma β (bottom two rows) in a top-down view (first and third rows) and an edge-on view(second and fourth rows). The white curves in the first two rows (panels a–j) are magnetic field streamlines, whereas in the last two rows (panels k–t) they are velocity vector field streamlines. Each column represents a different time, where t = 0 corresponds to the moment of protostellar birth. The mass of the protostar and its circumstellar disk are displayed in the top right corner of panels a–e. The scale bars in the first column apply to all the other columns as well.

Fig. 8. Same as Fig. 7, but for run IMHD. Since there is virtually no rotation in this run, the slices are done along the z-axis only. The green contour designates the stellar surface.

The magnetic field streamlines of run NIMHD are displayed in panels (a–j) of Fig. 7. In the top-down view (Fig. 7, panels a–e), we see the magnetic field lines being dragged by the nascent protostar. At the temperature-density regime displayed here, the ideal MHD limit is recovered as all dust grains are sublimated and ambipolar diffusion is no longer at play. In addition to this, the gas begins to ionize, which further enhances the magnetic field’s coupling to it. In addition, the plasma β values displayed in the last two rows indicates that thermal pressure support far outweighs magnetic pressure support, meaning that it is the fluid that dictates the behavior of the magnetic field. As such, the rotation of the newly-formed protostar and circumstellar disk causes a significant build-up of the toroidal component of the magnetic field, as the lines are twisted and tangled to an extreme degree by the dynamical second collapse. Along the disk midplane, there also seems to be a significant amount of turbulent magnetic eddies, which appear most prominent at later times. These are likely formed as a result of the emergence of spiral waves (see Figs. 3 and 4), which create significant turbulent motion within the disk. In essence, the turbulent eddies show that the magnetic field is at places confined within a tube-like structure which crosses the disk midplane, and hence showcases a significant poloidal component within the disk. We also note the spiral structure of the magnetic field intensity within the star-disk system. In Wurster et al. (2022), it is claimed that the Hall effect is responsible for the creation of this spiral structure; however, our results show here that the Hall effect is not necessary to form it. We show in Appendix D that the circumstellar disk is marginally stable against gravitational collapse, with the classical Toomre Q parameter hovering around unity.

Interestingly, in the edge-on view of panels (h–j) of Fig. 7, we see what appears to be a dipolar field in the western half of the star-disk system, where magnetic field streamlines originating from the southern pole of the protostar loop back into its northern pole; however, we are unsure as to why the feature appears only in the western half⁸. Outside the star-disk system, the magnetic field lines are mostly vertical and they thread the two bodies, showcasing the poloidal nature of the magnetic field in these regions. The plasma β decreases in the polar regions over time due to the depletion of material in these regions as the second collapse proceeds (Ahmad et al. 2024), which in turn causes a reduction in thermal pressure support. The disk’s surface also appears to have a plasma β ≈ 1, and the velocity vector field streamlines indicate that material is advected toward the protostar from the upper layers of the disk, as reported previously in the MHD run of Lee et al. (2021) and the hydro runs of Ahmad et al. (2024). Appendix E presents a quantitative measurement of the directional mass flux, which shows that the upper layers of the disk transport a similar amount of material toward the interior as the main body of the disk. Despite this, we see no outflow or high-velocity jet developing, as the velocity vector field streamlines in panels (p–t) of Fig. 7 are pointing toward the protostar, thus indicating infall. Any such outflows are likely to occur at much later times, when the polar reservoir of gas is significantly depleted and the plasma β in these regions drop to very small values. This once again confirms the results of Wurster & Lewis (2020b), which found that turbulence in the initial dense cloud core significantly delays the onset of jets and outflows. This is likely due to the absence of coherent magnetic field lines, which significantly hinders the onset of jets and outflows. In addition, Vaytet et al. (2018) also reports the absence of jets or outflows at protostellar scales despite the absence of initial turbulence in their progenitor dense core, which is likely due to the fact that the toroidal component of the magnetic field has yet to reach the strength needed to trigger the magneto-centrifugal mechanism (Lynden-Bell 1996; Lovelace et al. 2002). Finally, Machida et al. (2014) has shown a sensitivity of the presence of jets and outflows with regard to the large scale initial conditions where they struggled to recover the latter when using uniform density progenitor cores (as is used in Vaytet et al. (2018) and the present study).

With regard to the spatial distribution of the magnetic field within the star-disk system, we unsurprisingly see that the central region containing the protostar has the strongest field strength, reaching ≈5 × 10³ G. In accordance with Wurster et al. (2022)’s higher resolution runs, we witness spiral structures in magnetic field strength throughout the star-disk system.

In panels (h) and (i) of Fig. 6, quantitative measurements of B_ϕ and B_z are provided. The cylindrical radial velocity displayed in Fig. 6d allows one to locate the accretion shock, which manifests itself as a strong discontinuity⁹. First, at t = 0, the toroidal component is the dominant one within the protostar. However, at t ≈ 191 days, the vertical component is significantly built up and it becomes stronger than its toroidal counterpart. However, this appears to be transient, as the vertical component within the protostar seems to be oscillating. At larger radii (i.e., within the circumstellar disk), the opposite occurs: we see a build-up of the toroidal component of the magnetic field whereas the poloidal component is significantly reduced. In Fig. 6f, the black lines display the enclosed magnetic energy within the (spherical) radius r, which is computed as

$$\begin{array}{c}\hfill E_{B_{\mathrm{enc}}}=\frac{1}{2}{\displaystyle {\int }_0^r}{B}^2{r}^{2}dr.\end{array}$$(2)

We see that the innermost regions of the system lose magnetic energy over time. In these regions, the gas recovers the ideal MHD limit and flux freezing holds, with B ∝ ρ^2/3. Since the density within these regions remains somewhat constant, their loss of magnetic energy is due to an outward advection of material, as the protostar exceeds breakup velocity and begins shedding its surface material (Ahmad et al. 2024).

We now turn to describing the magnetic field structure of run IMHD (Fig. 8). Here, at t = 0 (panels (a) and (f) of Fig. 8), we see an extreme pinching of the magnetic field lines as a result of the second gravitational collapse. In addition, the field lines outside the protostar (lime contour) are almost entirely radial, with virtually no toroidal component present. However, as in Ahmad et al. (2023), we find strong turbulent motion within the protostar¹⁰. This causes all magnetic field components within the protostar to reach a similar strength, a fact that is particularly evident in panels (h) and (i) of Fig. 5.

One final result we would like to report is in regard to the slices shown in panels (f–j) of Fig. 8, displaying the plasma β of the gas. In Vaytet et al. (2018), it is reported that the protostar formed under the ideal MHD approximation is a magnetically supported object. However, we show here that all gas downstream of the protostellar accretion shock (lime contour) has a plasma β ≈ 1 or ≫1. This means that the protostar is an entirely thermally supported body¹¹. In addition, we have overlayed on these slices the velocity vector field streamlines, which show that no outflow or jet is being launched by the protostar. This is to be expected given the immensely unstructured nature of the magnetic field in this run, which exhibits no coherent toroidal component owing to the lack of rotational motion, and can thus no longer drive an outflow through the magneto-centrifugal mechanism (Blandford & Payne 1982; Ouyed & Pudritz 1997). This is in agreement with Wurster & Lewis (2020b).

3.4. Disk expansion: comparison with the hydro case

We now turn to providing a quantitative description of the evolution of the circumstellar disk over time. Since its structure appears to be qualitatively similar to the hydro case (a large and highly flared disk), we compare it to that obtained in run G2 of Ahmad et al. (2024), whose initial conditions and numerical setup are the same (notwithstanding the absence of magnetic fields). In addition, since Fig. 7 has shown that the plasma β ≫ 1 within the disk, we expect a similar evolution to the hydro case but with a notable increase in torquing owing to the strong magnetic field strength within the disk (∼[10² − 10³] G) and in the envelope (∼10² G). To this end, we display in Fig. 9 the mass, radius, specific angular momentum, and density at the equatorial shock front of the circumstellar disk in run NIMHD (black curves) and in the hydro run of Ahmad et al. (2024) (orange curves, hereafter run HD). We also leverage the information provided in Fig. 10, which displays the column density maps of runs NIMHD and HD (resp. panels (a) and (b) of Fig. 10) and their corresponding radial profiles (Fig. 10c). The velocity profiles are shown in Fig. 10d.

Fig. 9. Temporal evolution of the circumstellar disk of run NIMHD (black curves), compared with its hydro counterpart (orange curves; from Ahmad et al. 2024). The green curves are a zoom-out, branched from run NIMHD and run at a lower resolution with ℓmax = 24 (see Sect. 2.2). Panels (a), (b), (c), and (d) respectively display as a function of the disk’s equatorial radius Rd the density measured at the disk’s equatorial shock front (obtained through ray-tracing), the mass of the disk, the mass of the protostar, and the disk’s specific angular momentum. Panel (e) displays Rd as a function of time, where t = 0 marks the moment of birth of the disk. Panel (f) displays the specific angular momentum of the gas within 1 AU found inside the disk (solid and dashed lines) and outside the disk (dotted lines) as a function of time since the birth of the disk.

Fig. 10. Structural and kinematic comparison between run NIMHD and HD (resp. black and orange curves in panels c and d). Panel (a) displays column density maps for run NIMHD at our final simulation snapshot, which is ≈190 days after protostellar birth. Panel (b) displays the equivalent map for run HD, at a moment in time where its radius is comparable to that of run NIMHD (≈0.2 AU). Only cells belonging to the disk were used in the making of these maps. The second row displays the radial profiles of column density (panel c) and azimuthal veocity (panel d). The dashed lines in panel (d) display the Keplerian velocity computed with the protostar’s mass.

Additionally, since it is of interest to advance the simulation in time, we have branched run NIMHD at ≈0.4 years following protostellar birth and run a parallel simulation with a reduced ℓ_max = 24, which significantly alleviates the computational cost of the simulation. The properties of the disk in this run, labeled “NIMHD_LR”, are shown in the green curves of Fig. 9. The overlap between the black and green curves shows that its results are realistic enough for physical interpretations.

We first begin by studying the temporal evolution of the disk’s radius with respect to time (Fig. 9e). We note the fact that although the evolution in the initial 0.3 years are identical in the HD and NIMHD runs, they later diverge as run NIMHD exhibits a slower disk growth in time owing to strong torque mechanisms. Nevertheless, the plot displaying disk’s mass with respect to its radius (Fig. 9b) shows that the HD and NIMHD runs exhibit a similar evolutionary trend, although run NIMHD’s curves appears more oscillatory than that of run HD. These oscillations are caused by strong spiral waves within the disk in run NIMHD, which carry a significant amount of material inward and in doing so, also reduce the mass of the disk when compared to run HD. At a given disk radius, run NIMHD displays a smaller mass than run HD. These spiral waves are likely caused by magnetic torques, which reduce the gas’s centrifugal support against gravity, and they warp the disk (Lai 2003; Tomida et al. 2013), as can be seen in Fig. 10a. When measuring the mass-weighted mean magnetic field strength within a radius of 10 AU, the toroidal component outweighs all others by an order of magnitude (≈270 G, compared to ≈60 G for the cylindrical radial component and ≈25 G for the vertical field). This means that the magnetic field is mostly parallel to the disk midplane, which has been shown to cause prominent spiral waves to develop (Joos et al. 2012; Li et al. 2013; Hennebelle et al. 2020b). In contrast, the hydro disk (Fig. 10b) remains circular. This causes higher column densities in the outer regions of the disk in run HD than in run NIMHD (Fig. 10c). Finally, the higher protostellar mass causes faster rotation in the innermost regions of the disk in run NIMHD (Fig. 10d), and its velocity profile closely approaches the Keplerian profile (dashed line), whereas in run HD the velocity profile is super-Keplerian owing to the disk’s self-gravity. The disk is slightly sub-Keplerian in run NIMHD owing to thermal pressure support. Magnetic pressure support has a negligible effect on v_ϕ since β ≪ 1.

The prominence of magnetic braking is further displayed in Fig. 9c, which shows the protostellar mass as a function of disk radius. Here, we see that although run HD quickly reaches a plateau in disk mass, run NIMHD shows a rapidly growing protostar. We note that the specific angular momentum of the disk, shown in Fig. 9d, is the same in both runs and scales as $\sqrt{R_{\mathrm{d}}}$¹². In Fig. 9f, we show the specific angular momentum of the gas within 1 AU, computed both outside of the disk (dotted lines) and within it (solid and dashed lines). This figure shows that the disk in run NIMHD is accreting from a reservoir of gas containing a smaller amount of angular momentum than run HD, which shows that magnetic braking occurred before the gas was accreted onto the disk.

The reduced mass of the disk in run NIMHD when compared to run HD also manifests itself in a reduced disk density. More specifically, when measuring the density at the disk’s equatorial shock front (ρ_acc, Fig. 9a), we see that it is consistently lower than in run HD. Although we could not integrate the calculations to the point where the disk reaches the commonly used sink accretion radius of 1 AU, we present in Appendix C the results of a run branched out of NIMHD_LR (labeled NIMHD_LR_2), whose lower resolution allowed for the calculation to be pushed longer out in time. It predicts ρ_acc(1 AU)≈2 × 10⁻¹⁰ g cm⁻³, a value that is a factor if ≈3 lower than that predicted by the hydrodynamical run (≈5.93 × 10⁻¹⁰ g cm⁻³).

3.5. Eccentricity measurements

Recently, Commerçon et al. (2024) has shown that disks born out of gravitational collapses exhibit strong eccentricity. In the present study, run NIMHD has yielded a disk that also appears to show eccentric motion, as is most prominently seen in the asymmetrical flaring of the disk in panels (p–t) of Fig. 7. As such, we provide quantitative measurements of the eccentricity of the circumstellar disk in run NIMHD. In order to do so, we compute both the eccentricity vector e and semimajor axis a, defined as (Commerçon et al. 2024)

$$\begin{array}{cc}\hfill \mathit{e}& =\frac{\mathit{v}\times \mathit{j}}{G{M}_{\mathrm{enc}}(r)}-\frac{\mathit{r}}{r},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill a& =\frac{rG{M}_{\mathrm{enc}}(r)}{2G{M}_{\mathrm{enc}}(r)-v^{2}r},\hfill \end{array}$$(4)

where j is the specific angular momentum of the gas and r is the position vector. Since the disk mass is greater than or comparable to the protostellar mass (see panels (b) and (c) of Fig. 9), we compute e and a using M_enc(r), the enclosed mass within r:

$$\begin{array}{c}\hfill M_{\mathrm{enc}}=4\pi {\displaystyle {\int }_0^r}\rho r^{2}dr.\end{array}$$(5)

The resulting measurements of the magnitude of the eccentricity vector e = |e| are presented in Fig. 11, which displays the mass-weighted orbital averages of e. We see in this figure strong eccentric motions within the disk, with e reaching values that consistently exceed 0.2. This shows that similarly to larger scale disks forming around sink particles (e.g., as in Lee et al. 2021; Commerçon et al. 2024), the circumstellar disk forming as a result of the rotational breakup of the protostar is also eccentric. The strong gradient in e in the outer regions of the disk, seen particularly at later times, are caused by strong spiral motions within the disk.

Fig. 11. Mass-weighted average of the magnitude of the eccentricity vector as a function of semimajor axis (see Eqs. (3) and (4)) in run NIMHD. Each curve corresponds to a different time, where t = 0 corresponds to the epoch of protostellar birth.

4. Discussions

The results presented here, most notably those of run NIMHD, are noteworthy for a number of outstanding issues in stellar formation theory. We discuss their implications in this wider context below.

4.1. The angular momentum problem

In Ahmad et al. (2024), we reported on the birth of the circumstellar disk as a result of the breakup of the protostar and the subsequent vigorous radial expansion of the disk in time. Run NIMHD has confirmed that such a phenomenon is reproduced even in the presence of magnetic fields, provided that magnetic resistivities are accounted for. In the literature, Machida & Matsumoto (2011), Vaytet et al. (2018), Wurster & Lewis (2020b), Machida & Basu (2019), Bhandare et al. (2024) have also reported on the birth of a circumstellar disk that rapidly expands to larger radii. Machida et al. (2007) also reported that their non-ideal MHD calculations produced rapidly rotating protostars, and emphasized that the calculations must be run farther out in time to properly model the angular momentum evolution of prestellar objects.

What these simulations seem to show is that a paradigm shift is required in our understanding of the angular momentum problem. It has long been implicitly implied in stellar formation theory that angular momentum must be lost during the collapse so as to prevent the protostar from ever reaching breakup velocity (Bodenheimer 1995). What our comparison between run IMHD and NIMHD shows is that should angular momentum transport by magnetic fields be so efficient so as to prevent the protostar from ever reaching breakup velocity, then no circumstellar disk forms: it is the very fact that parcels of fluid within the protostar achieve rotational breakup that allows for circumstellar disks to form in our simulations¹³. As such, whatever angular momentum transport process is responsible for spinning down the protostar to the ∼10 − 15% of breakup velocity as is observed in YSOs (Rebull et al. 2004; Herbst et al. 2007), it is acting on longer timescales than the free-fall time of the cloud.

The vigorous radial expansion of the disk may even have left its mark on the meteoric record of the solar system, as high-temperature condensates in the form of calcium-aluminum inclusions and amoeboid olivine aggregates show evidence of rapid outward transport (Morbidelli et al. 2024).

4.2. The magnetic flux problem

As mentioned previously, YSOs are consistently found to have ∼kG magnetic field strengths, with the earliest measurements being those of Class I sources. It is as of yet unclear what the origin of these magnetic fields is; whether they are fossil fields that are maintained for ∼10⁵ yr, or generated during the pre-main sequence evolutionary phase through a dynamo process. Our results, and those of Vaytet et al. (2018), Machida & Basu (2019), Wurster et al. (2022)¹⁴, show that the measured values may be achieved and maintained following the second collapse. The uncertainty here lies in the short horizon of predictability of second collapse calculations, as they must be able to simulate ∼10⁵ years following protostellar birth in order to make an adequate comparison with observations. As such, in order to determine the origin of the magnetic fields measured in YSOs, better constraints on magnetic field strengths both at dense core scales through the measurement of linearly polarized dust emissions, or at much smaller scales through Zeeman broadening, are required in Class 0 sources. However, these measurements are immensely difficult to undertake owing to the optical depths involved.

In the meantime, a significant amount of theoretical modeling is required in order to describe the evolution of the protostellar magnetic field in conjunction with pre-stellar evolution models, however simplified such models are. Numerically costly simulations such as those presented in this study are immensely helpful in obtaining the initial properties and structure of the protostar; however, their short horizon of predictability precludes them from definitively solving the magnetic flux problem.

4.3. The missing mass problem

Current observational surveys of Class 0/I disks estimate their masses to be ∼10⁻³ − 10⁻² M_⊙ (e.g., Tobin et al. 2020), which appears to be an order of magnitude lower than those predicted by theoretical studies (∼10⁻² − 10⁻¹ M_⊙; e.g., Machida & Matsumoto 2011; Tsukamoto et al. 2015b,a; Tomida et al. 2015; Masson et al. 2016; Lee et al. 2021, see the discussions in Tsukamoto et al. 2023). Notwithstanding the uncertainties involved in current observational methods (Tung et al. 2024), this discrepancy has been dubbed the “missing mass problem”. It has also been shown that current subgrid models aiming to emulate the sub-AU regions by replacing them with a sink particle show a strong sensitivity to the parameters chosen in this model. In Hennebelle et al. (2020a), the sink accretion threshold was shown to particularly affect the disk mass, with lower accretion thresholds leading to lower disk masses. The results of the hydrodynamical runs of Ahmad et al. (2024) seemed to show that the sink accretion threshold used in most simulations (≈1.66 × 10⁻¹¹ g cm⁻³) is a factor of ≈40 lower than it should be, thus exacerbating the missing mass problem as that would mean that disks are in reality much more massive in simulations. In the present study, the inclusion of magnetic fields results in a lower disk density, and by extension a lower density measured at the disk shock front, placing it at ≈2 × 10⁻¹⁰ g cm⁻³ when it reaches a radius of 1 AU, as opposed to the ≈5.93 × 10⁻¹⁰ g cm⁻³ obtained in the hydrodynamical run. This result highlights the sensitivity of the nascent circumstellar disk to the magnetic field strength that is inherited by the sub-AU region. Thus, a more thorough understanding of the disk properties not only requires one to study the nascent protostar and disk in concert, but also requires better constraints on dust resistivities that dictate the amount of magnetic flux, and by extension the amount of angular momentum inherited by the disk.

In any case, this issue highlights the need for better comparisons between observations and theoretical models, as Tung et al. (2024) has shown that current observational estimates of disk masses are inadequate and fail to predict the current sizes when compared to simulations. Advancements in this regard are of great importance to the field, as constraining the masses of circumstellar disks is crucial to determine the initial mass budget for planet formation.

4.4. The importance of adequate dust resistivity tables

An important uncertainty in our current understanding of sub-AU regions is the dust resistivity used. The MRN dust size distribution is increasingly called into question by studies that account for dust coagulation and fragmentation during protostellar collapses (Lebreuilly et al. 2023; Kawasaki & Machida 2023; Tsukamoto et al. 2023; Bhandare et al. 2024; Vallucci-Goy et al. 2024). Our simulations are undertaken under the assumption that Ohmic resistivity is, as predicted by dust-size distribution studies, negligible. This leads to stronger magnetic fields within the first Larson core, which in turn increases the magnetic field intensity in the nascent protostar and circumstellar disk. As a result, magnetic torques drive considerably more material toward the protostar, thus leading to a reduced disk density.

As such, the properties of the newly-formed circumstellar disk are highly sensitive to the dust resistivities that dictate the magnetic field intensity inherited from larger spatial scales. A better understanding of the sub-AU regions is predicated upon accurate dust resistivity tables, which requires a better understanding of the dust size distribution. The Hall effect, however, remains an important caveat as there appears to be conflicting models in the literature regarding its effects during the collapse (Wurster 2021; Tsukamoto & Okuzumi 2022; Hopkins et al. 2024). Progress will ultimately be achieved by longer wavelength observations of star-forming regions in order to probe optically thin dust emissions, as well as by advances in our theoretical modeling of dust growth and fragmentation during protostellar collapses. Recent numerical advances in the description of these processes may allow for their inclusion in fully 3D hydrodynamical simulations (Lombart et al. 2024).

4.5. A magneto-rotational instability?

In run NIMHD, we witness a hot and highly magnetized disk (T ∼ 10³ K and B ∼ [10² − 10³] G), see Figs. 3, 6 and 7), whose plasma β, as shown in Fig. 7 is well above unity. This naturally calls for an analysis to determine whether such a disk is prone to develop the magneto-rotational instability (MRI, Balbus & Hawley 1991; Wardle 1999; Lesur et al. 2014; Riols & Latter 2019; Deng et al. 2020), an instability known to occur in such conditions. Although we chose to omit Ohmic dissipation from our simulation, the analysis hereafter presented also accounts for Ohmic resistivity, as it is the non-ideal MHD effect commonly responsible for quenching MRI in other studies in the literature (Lesur et al. 2014). Hence, its inclusion in this section’s analysis is purely to determine whether or not it may quench the MRI were we to include it.

In order to ascertain whether an MRI can operate in run NIMHD, we first evaluate whether the ideal MHD limit is recovered within the star-disk system. This can be done by computing the ambipolar and Ohmic Elsässer numbers (respectively Λ_AD and Λ_O)

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_{\mathrm{AD}}& =\frac{v_{\mathrm{A},\mathrm{z}}^2}{{\eta }_{\mathrm{AD}}\mathrm{\Omega }},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\mathrm{\Lambda }}_{\mathrm{O}}& =\frac{v_{\mathrm{A},\mathrm{z}}^2}{{\eta }_{\mathrm{O}}\mathrm{\Omega }},\hfill \end{array}$$(7)

where v_A, z is the vertical Alfvén velocity, η_AD (respectively η_O) is the ambipolar (respectively Ohmic) resistivity¹⁵, and Ω is the angular velocity. When the Elsässer number is maintained below unity, the MRI is effectively quenched (Jin 1996; Turner et al. 2007; Simon et al. 2013).

The resulting measurements of Λ_AD and Λ_O are presented in Fig. 12, which displays the two quantities as a function of radius within the protostar and circumstellar disk in 2D histograms. The intensity of the color shade is linked to the average density within each bin. As a complement to this figure, we show in Fig. 13 both top-down (frist row) and edge-on (second row) slices displaying Λ_AD (first column) and Λ_O (second column).

Fig. 12. Set of 2D histograms binning the cells belonging to the protostar and circumstellar disk (see Appendix A for an overview of each object’s definition) at the final simulation snapshot of run NIMHD (t ≈ 0.55 yr after protostellar birth). These display the ambipolar (ΛAD, blue distribution) and Ohmic (ΛO, green distribution) Elsässer numbers as a function of radius. The gray dotted line denotes Λ = 1. The color shading is related to the average density within the bin, as indicated by the grayscale colorbar.

Fig. 13. Slices showing the ambipolar (first column, panels a and c) and Ohmic (second column, panels b and d) Elsässer numbers within the star-disk system in the final snapshot of run NIMHD (t ≈ 0.55 yr after protostellar birth). The first row (panels a and b) displays slices in a top-down view, whereas the second row (panels c and d) displays edge-on views. The colorbar on the left applies to panels (a) and (c); the colorbar on the right applies to panels (b) and (d). The outer contour corresponds to the disk surface, whereas the inner one corresponds to the protostar’s surface (see Appendix A for an overview of each object’s definition).

We see that Λ_AD is above unity throughout the system, reaching values of ∼10¹⁶. It is below unity only for low-density gas at larger radii, which appear to be rather sparsely distributed, as seen in panels (a) and (c) of Fig. 13. The value of Λ_O, however, is consistently above unity only in the innermost regions belonging to the protostar, and it fluctuates throughout the rest of the star-disk system according to the local magnetic field strength¹⁶. These measurements show that the ideal MHD limit is mostly recovered within the star-disk system, and should we have included Ohmic dissipation in our calculations, it would have been mostly recovered in the innermost regions of the system.

In order to assess whether the conditions to trigger the MRI are met, one final criterion is necessary, a comparison of the disk scale height $h=\frac{c_{\mathrm{s}}}{\mathrm{\Omega }}$ (where c_s is the sound speed) to the wavelength of the most unstable MRI mode (λ_MRI, Balbus & Hawley 1991):

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{MRI}}=\frac{2\pi }{\sqrt{3}}\frac{v_{\mathrm{A},\mathrm{z}}}{\mathrm{\Omega }}·\end{array}$$(8)

The MRI operates when h > λ_MRI (Balbus & Hawley 1991), a condition that is met within the disk as h ∈ [10⁻³, 10⁻¹] AU, whereas λ_MRI ∈ [10⁻¹⁰, 10⁻²] AU. As such, we conclude that the circumstellar disk born out of the second gravitational collapse is prone to triggering a dynamo process through the MRI mechanism.

However, given the small values of λ_MRI witnessed within the star-disk system, we do not have enough resolution to adequately resolve the instability. This is shown in Fig. 14, which displays λ_MRI/Δx in slices at the final snapshot of run NIMHD. λ_MRI/Δx is consistently below unity throughout most of the circumstellar disk, and it fails to reach 10 in the midplane, which is insufficient to adequatly capture the instability (Noble et al. 2010). The only regions where λ_MRI/Δx ≫ 1 are the polar regions directly above the protostar, where the angular velocity is smallest. This could be behind the previously discussed oscillations in magnetic field strength in the innermost regions.

Fig. 14. Top-down (panel a) and edge-on (panel b) slices across the star-disk system showing the ratio of the wavelength of the most unstable MRI mode to cell size at the final snapshot of run NIMHD (t ≈ 0.55 yr after protostellar birth). The outer contour corresponds to the disk surface, whereas the inner one corresponds to the protostar’s surface.

In the present study, we have shown strong torque mechanisms on the circumstellar disk that are driving high mass-accretion rates toward the protostar, even without the presence of the MRI. Nevertheless, the fact that the MRI may operate within this circumstellar disk is a very significant result, as it may significantly alter the structure of the magnetic field, as well as amplify its strength to equipartion values. This could later drive a high-velocity jet from protostellar scales, as well as induce even stronger torques on the disk, thus causing a higher accretion rate toward the protostar. This might ultimately cause the protostar to decouple from the circumstellar disk, and truncate the latter at a magnetospheric radius, as is observed in more evolved class I systems. Capturing the MRI in a circumstellar disk formed out of the second collapse might not be viable with current computational hardware, as the characteristic scales are too difficult to resolve. However, a different numerical setup inspired by the results obtained in this study might achieve such a result.

5. Conclusion

In the present study we presented radiative MHD simulations describing the collapse of a turbulent and gravitationally unstable dense cloud core of 1 M_⊙ to stellar densities, under the ideal MHD approximation and under the non-ideal approximation in which we accounted for ambipolar diffusion. Our stringent refinement criterion, as well as high spatial resolution, allowed us to describe the nascent protostar and circumstellar disk with unprecedented resolution. We pushed the calculations as far as possible in time following protostellar birth in order to study the nascent disk’s expansion, reaching ≈0.5 years in our high-resolution run and ≈1.2 years in our lower resolution run. Our results may be summarized as follows:

• (i)

  When accounting for ambipolar diffusion, the efficacy of magnetic braking is significantly reduced toward higher density gas, which allows the nascent protostar to reach breakup velocity and shed its surface material to form a circumstellar disk around it. The nascent disk exhibits strong eccentricity (reaching values of ∼0.1), and the protostar is embedded within it. The birth and early evolution of the circumstellar disk is qualitatively similar to the RHD runs presented in Ahmad et al. (2024), as the plasma β within the disk far exceeds unity. The nascent disk vigorously expands in the radial direction. This result carries implications for the angular momentum problem, as we show that the protostar must achieve breakup velocity in order to form a circumstellar disk. As such, angular momentum transport processes must spin-down the protostar on considerably longer timescales than the free-fall time of the dense cloud core.

• (ii)

  The magnetic field implanted in the protostar at birth has a strength of ∼10⁵ G in the ideal MHD run, which then continuously reduces to ∼10⁴ G as the simulation progresses. In the non-ideal MHD run, the implanted field has a strength of ∼10³ G which is maintained throughout the simulation’s duration. Since current observational surveys of magnetic fields in YSOs favor the fossil field hypothesis, this puts the non-ideal MHD simulation in agreement with them.

• (iii)

  The field implanted in the protostar in the non-ideal run is mostly toroidal (B_ϕ), although a notable vertical component (B_z) threads the star-disk system. Within the protostar, the vertical component is significantly built up over time.

• (iv)

  The circumstellar disk formed in our non-ideal run has a plasma β well above unity, with a strong magnetic field whose strength is in the range [10² − 10³] G. Our analysis shows that this disk is prone to triggering a dynamo process through the magneto-rotational instability (MRI), although we do not have the resolution to adequately capture the mechanism. This would be the case even if we were to account for Ohmic dissipation. The MRI might be responsible for the decoupling of the protostar from the disk in which it is embedded, and transition the system toward a magnetospheric accretion mechanism reminiscent of class I systems.

• (v)

  Owing to our use of turbulent initial conditions, the magnetic field mostly loses its coherence and we see no outflows or jets in both runs. In the non-ideal MHD run, however, the plasma β in the polar regions upstream of the protostellar accretion shock is continuously being reduced. Coupled with the fact that a vertical component is being built up in the protostar, this may lead to the launching of an outflow at later times.

• (vi)

  When comparing the nascent disk in the non-ideal MHD run to its hydro counterpart (from Ahmad et al. 2024), we note a reduced disk density. This is caused by the presence of strong magnetic and gravitational torques within the disk that transport a significant amount of material toward the protostar. This also causes the protostar to become more massive than in the hydro case. The reduced disk density in turn causes a reduced density at the disk’s equatorial shock front, which is an important measure for studies of global disk evolution that leverage sink particles to advance the simulations in time. The trends seen in our simulations indicate that the shock front’s density at 1 AU in the magnetized case is a factor of ≈3 lower than that reported in the hydro runs of Ahmad et al. (2024). As such, we conclude that constraining the current subgrid model parameters used in the literature requires better constraints on dust resistivities, thus highlighting the need for a more comprehensive modeling of the dust size distribution throughout the collapse.

Although we may learn a great deal from expensive simulations like the ones presented in this study, it is important to note that their horizon of predictability is rather short and that their results may not be applicable throughout the entirety of the Class 0 phase. The importance of magnetic fields in dictating the transport of material within the circumstellar disk also highlights the need for better constraints on the dust-size distribution, which requires significant observational and theoretical efforts.

Acknowledgments

We thank the anonymous referee for their numerous comments that have significantly improved the quality of the manuscript. This work has received funding from the French Agence Nationale de la Recherche (ANR) through the projects COSMHIC (ANR-20-CE31-0009), DISKBUILD (ANR-20-CE49-0006), and PROMETHEE (ANR-22-CE31-0020). We have also received funding from the European Research Council synergy grant ECOGAL (Grant : 855130). The simulations were carried out on the Alfven super-computing cluster of the Commissariat à l’Énergie Atomique et aux énergies alternatives (CEA). We thank Elliot Lynch for valuable discussions during the writing of this paper. Post-processing and data visualization was done using the open source Osyris package.

References

Appendix A: Defining the protostar and circumstellar disk

Since the two runs yield drastically different qualitative properties, one must have robust definitions for both the protostar and the circumstellar disk before proceeding to any quantitative comparison. In the case of run IMHD, the protostar is a thermally supported body; however, simply finding the cells in which thermal pressure support is attained as is done in Ahmad et al. (2023) is inadequate, as the current sheets protruding from the stellar surface also satisfy this definition. As such, we have decided to define the protostar as being all cells in which at least 90% of H₂ molecules are dissociated (i.e., X_H2 < 10⁻¹, where X_H2 is the fraction of hydrogen under molecular form).

In the case of run NIMHD, the presence of centrifugal support drastically changes the structure of the protostar, which flattens along the equator. Further complicating things, the transition from thermal pressure support to mainly centrifugal support against gravity is smooth, and no shock front separates the protostar from its circumstellar disk (Ahmad et al. 2024). As a result, we adopt the same arbitrary definition for the protostar as in Vaytet et al. (2018), Ahmad et al. (2024), namely, that it is the gas whose density exceeds 10⁻⁵ g cm⁻³. To illustrate why these two criteria were used, we display in Fig. A.1 their results when applied to both simulations. The criterion defining the protostar as being X_H2 < 10⁻¹ is displayed in the first row, where we see that it recovers the stellar surface in run IMHD but fails to do so in run NIMHD. On the other hand, the second criterion stating that the star is defined as ρ > 10⁻⁵ g cm⁻³ and displayed in the second row, shows that it selects extended current sheets protruding from the stellar surface in run IMHD but recovers a centrifugally flattened surface in run NIMHD.

Fig. A.1. 3D illustration of the two criteria used to define the protostar, as applied in run IMHD (first column) and run NIMHD (second column). The first row displays an isocontour of XH2 ≈ 10−1, whereas the second row displays an isocontour of ρ ≈ 10−5 g cm−3. The colorbar in the first (resp. second) row displays the gas density (resp. XH2) in the extracted surface.

The circumstellar disk is defined as in Ahmad et al. (2024); it is the centrifugally supported gas whose thermal pressure support exceeds incoming ram pressure, and whose density exceeds the density of the shock front (which is in turn determined through ray-tracing).

Appendix B: Rotational breakup of the protostar

Here we present the measurements of $v_{\varphi }/\sqrt{|g_{\mathrm{r}}|r}$, where g_r = −∂ϕ/∂r (ϕ being the gravitational potential obtained through the Poisson equation), as is done in Ahmad et al. (2024), in order to demonstrate that the protostar in run NIMHD undergoes a rotational breakup. The results, presented in Fig. B1, show that parcels of fluid exceed breakup velocity, thus causing the gas to spread outward due to excess angular momentum.

Fig. B1. 2D histogram binning all the cells in run NIMHD at t ≈ 0.3 yr (where t = 0 corresponds to the epoch of protostellar birth), which shows the distribution of azimuthal velocities divided by $\sqrt{|g_{\mathrm{r}}|r}$ with respect to radius. The color-coding in the histograms represents the co-latitude θ divided by π, where θ/π = 0.5 corresponds to the equator and θ/π = 1 (respectively 0) corresponds to the south (respectively north) pole.

Appendix C: An additional zoom-out to predict ρ_acc(1 AU)

One of the goals of this study is to provide a prediction for ρ_acc(1 AU), the disk’s equatorial shock front density when it reaches a size of 1 AU. Run NIMHD_LR, although pushed much farther out in time than run NIMHD, is still prohibitively expensive. As such, we have branched run NIMHD_LR at t ≈ 0.92 yr after protostellar birth and pushed the calculations until the disk radius reached ≈1 AU. The results of this run, labeled NIMHD_LR_2 are displayed in the red curve of Fig. C.1.

Fig. C1. Density measured at the disk’s equatorial shock front as a function of the disk’s equatorial radius Rd, for run NIMHD (black curve) and its hydro counterpart (orange curve). The green curve is a zoom-out, branched from run NIMHD and run at a lower resolution with ℓmax = 24 (run NIMHD_LR), and the red curve is a zoom-out, branched out of run NIMHD_LR and run at ℓmax = 23 (run NIMHD_LR_2). The significant overlap between the results of NIMHD_LR and NIMHD_LR_2 (green and red curves) shows that the results of this doubly zoomed-out run are converged with regard to ρacc.

Appendix D: Gravitational stability of the circumstellar disk

In this section we provide measurements of the Toomre Q parameter (Toomre 1964)

$$\begin{array}{c}\hfill Q=\frac{\omega c_{\mathrm{s}}}{\pi G\mathrm{\Sigma }},\end{array}$$(D.1)

where ω is the epycyclic frequency. This parameter analyzes the gravitational stability of the circumstellar disk formed in run NIMHD, with Q > 1 indicating a disk that is stable against gravitational instabilities.

The resulting measurements are presented in Fig. D1, which shows that Q ∼ 1 throughout most of the circumstellar disk, thus indicating marginal stability, as is expected of self-gravitating accretion disks (Lodato & Rice 2004; Lodato 2007).

Fig. D1. Real part of the Toomre Q parameter of the circumstellar disk of run NIMHD, seen in a top-down view where each panel corresponds to a different time, with t = 0 corresponding to the epoch of protostellar birth. Only cells belonging to the circumstellar disk were selected when producing these plots.

Appendix E: Directional mass flux

In this section we provide measurements of the radial mass flux within the circumstellar disk in run NIMHD, both for the upper layers of the disk, and its main body. These are defined as:

$$\begin{array}{c}\hfill \{\begin{array}{cc}\text{Upper layers: Disk cells such that}\rho <3{\rho }_{\mathrm{S}},\hfill & \\ \text{Main body: All remaining disk cells},\hfill \end{array}\end{array}$$

where ρ_S is the disk’s shock density obtained through ray-tracing. This definition was chosen such that the 3ρ_S contour produces a visually compelling result for the upper layers of the disk, comprising of the region just downstream of the shock front. The radial mass flux is computed using −ρv_r, which is then averaged within each radial bin. The resulting measurements are presented in Fig. E1 massflux, which show that material is flowing inward throughout most the disk, both in the upper layers and in the main body. The innermost regions of the main body of the disk show outward transport due to excess angular momentum.

Fig. E1. Average radial mass flux within the circumstellar disk in run NIMHD, measured in radial bins at t ≈ 0.5 yr (where t = 0 corresponds to the epoch of protostellar birth). Measurements are performed along the upper layers of the disk (solid line) and the main body (dashed line). Only cells belonging to the circumstellar disk were used in the computation.

All Figures

Fig. 1. Comparison of runs IMHD (first row) and NIMHD (second row) at the scale of the dense cloud core itself. The snapshots are taken respectively at t ≈ 23.25 and t ≈ 23.39 kyr following the collapse of the dense core. The first column displays column density (panels a and b), the second column displays the optical depth computed along the line of sight (panels c and d), and the last column displays the maximum temperature along the line of sight (panels e and f). All maps are projections along the z-axis. The green contour in panels (c) and (d) represent an optical depth of unity. The scale bars in the first column apply to the other columns as well.

In the text

Fig. 2. Quantitative comparison of the collapse between run IMHD (black) and run NIMHD (red). Panel (a) displays the evolution of the maximum density since the formation of the first Larson core, which we define as the time when a density of ≈10−10 g cm−3 is achieved. Panels (b) and (c) display the magnetic field strength evolution as a function of time since first core formation and since protostellar birth (defined as the moment a density of ≈10−5 g cm−3 is reached), where the solid lines represent the maximum magnetic field strength and the dotted lines represent the field’s strength measured at the location of maximum density.

In the text

Fig. 3. Set of slices showing the evolution of the density (first two columns) and temperature (last two columns) for run IMHD (first and third columns, panels a–e and k–o) and run NIMHD (second and fourth columns, panels f–j and p–t). Each row represents a different time, where t = 0 corresponds to the moment of protostellar birth. For comparative purposes, the slices are shown at similar times, and the timestamp is written in the top right corner of each panel. The slices are done in the z direction for run IMHD, and along the angular momentum vector for run NIMHD. The scale bars in the first row apply to all the other rows as well.

In the text

	Fig. 4. Same as Fig. 3, but now showing the radiative flux (first two columns) and radial velocity (last two columns).
In the text

Fig. 5. Set of measurements of various physical properties of run IMHD at protostellar birth (t = 0, solid lines) and t ≈ 117 days later (dotted lines). These are averages in spherical bins, and show the gas density (panel a), temperature (panel b), entropy (panel c), radial and azimuthal velocity (panels d and e), magnetic field intensity (panel g), and the azimuthal and meridional components of the magnetic field, normalized by its magnitude (panels h and i). Panel (f) displays the enclosed magnetic energy (black lines) and the enclosed mass (red lines) as a function of spherical radius.

In the text

Fig. 6. Set of measurements of various physical properties of run NIMHD. These were done in cylindrical radial bins, in which only cells belonging to the midplane, defined as z ∈ [ − 2.5;  2.5]×10−2 AU, were used. The solid lines are measurements done at t = 0 (corresponding to the moment of protostellar birth), and the dotted lines are measurements done ≈190 days later. Panels (a), (b), and (c) display respectively the gas density, temperature, and specific entropy. Panels (d) and (e) display the gas (cylindrical) radial and azimuthal velocity. The red curves in panel (e) display the Keplerian velocity computed with the protostar’s mass at a given snapshot. Panels (g), (h), and (i) display the magnetic field strength, its toroidal component, and its vertical component. The toroidal and vertical components are normalized by the total magnetic field strength. Panel (f) displays the enclosed magnetic energy (black lines) and the enclosed mass (red lines) as a function of spherical radius.

In the text

Fig. 7. Magnetic field of run NIMHD: a set of slices displaying the magnetic field strength (top two rows) and plasma β (bottom two rows) in a top-down view (first and third rows) and an edge-on view(second and fourth rows). The white curves in the first two rows (panels a–j) are magnetic field streamlines, whereas in the last two rows (panels k–t) they are velocity vector field streamlines. Each column represents a different time, where t = 0 corresponds to the moment of protostellar birth. The mass of the protostar and its circumstellar disk are displayed in the top right corner of panels a–e. The scale bars in the first column apply to all the other columns as well.

In the text

	Fig. 8. Same as Fig. 7, but for run IMHD. Since there is virtually no rotation in this run, the slices are done along the z-axis only. The green contour designates the stellar surface.
In the text

Fig. 9. Temporal evolution of the circumstellar disk of run NIMHD (black curves), compared with its hydro counterpart (orange curves; from Ahmad et al. 2024). The green curves are a zoom-out, branched from run NIMHD and run at a lower resolution with ℓmax = 24 (see Sect. 2.2). Panels (a), (b), (c), and (d) respectively display as a function of the disk’s equatorial radius Rd the density measured at the disk’s equatorial shock front (obtained through ray-tracing), the mass of the disk, the mass of the protostar, and the disk’s specific angular momentum. Panel (e) displays Rd as a function of time, where t = 0 marks the moment of birth of the disk. Panel (f) displays the specific angular momentum of the gas within 1 AU found inside the disk (solid and dashed lines) and outside the disk (dotted lines) as a function of time since the birth of the disk.

In the text

Fig. 10. Structural and kinematic comparison between run NIMHD and HD (resp. black and orange curves in panels c and d). Panel (a) displays column density maps for run NIMHD at our final simulation snapshot, which is ≈190 days after protostellar birth. Panel (b) displays the equivalent map for run HD, at a moment in time where its radius is comparable to that of run NIMHD (≈0.2 AU). Only cells belonging to the disk were used in the making of these maps. The second row displays the radial profiles of column density (panel c) and azimuthal veocity (panel d). The dashed lines in panel (d) display the Keplerian velocity computed with the protostar’s mass.

In the text

	Fig. 11. Mass-weighted average of the magnitude of the eccentricity vector as a function of semimajor axis (see Eqs. (3) and (4)) in run NIMHD. Each curve corresponds to a different time, where t = 0 corresponds to the epoch of protostellar birth.
In the text

Fig. 12. Set of 2D histograms binning the cells belonging to the protostar and circumstellar disk (see Appendix A for an overview of each object’s definition) at the final simulation snapshot of run NIMHD (t ≈ 0.55 yr after protostellar birth). These display the ambipolar (ΛAD, blue distribution) and Ohmic (ΛO, green distribution) Elsässer numbers as a function of radius. The gray dotted line denotes Λ = 1. The color shading is related to the average density within the bin, as indicated by the grayscale colorbar.

In the text

Fig. 13. Slices showing the ambipolar (first column, panels a and c) and Ohmic (second column, panels b and d) Elsässer numbers within the star-disk system in the final snapshot of run NIMHD (t ≈ 0.55 yr after protostellar birth). The first row (panels a and b) displays slices in a top-down view, whereas the second row (panels c and d) displays edge-on views. The colorbar on the left applies to panels (a) and (c); the colorbar on the right applies to panels (b) and (d). The outer contour corresponds to the disk surface, whereas the inner one corresponds to the protostar’s surface (see Appendix A for an overview of each object’s definition).

In the text

	Fig. 14. Top-down (panel a) and edge-on (panel b) slices across the star-disk system showing the ratio of the wavelength of the most unstable MRI mode to cell size at the final snapshot of run NIMHD (t ≈ 0.55 yr after protostellar birth). The outer contour corresponds to the disk surface, whereas the inner one corresponds to the protostar’s surface.
In the text

	Fig. A.1. 3D illustration of the two criteria used to define the protostar, as applied in run IMHD (first column) and run NIMHD (second column). The first row displays an isocontour of XH2 ≈ 10−1, whereas the second row displays an isocontour of ρ ≈ 10−5 g cm−3. The colorbar in the first (resp. second) row displays the gas density (resp. XH2) in the extracted surface.
In the text

Fig. B1. 2D histogram binning all the cells in run NIMHD at t ≈ 0.3 yr (where t = 0 corresponds to the epoch of protostellar birth), which shows the distribution of azimuthal velocities divided by $\sqrt{|g_{\mathrm{r}}|r}$ with respect to radius. The color-coding in the histograms represents the co-latitude θ divided by π, where θ/π = 0.5 corresponds to the equator and θ/π = 1 (respectively 0) corresponds to the south (respectively north) pole.

In the text

Fig. C1. Density measured at the disk’s equatorial shock front as a function of the disk’s equatorial radius Rd, for run NIMHD (black curve) and its hydro counterpart (orange curve). The green curve is a zoom-out, branched from run NIMHD and run at a lower resolution with ℓmax = 24 (run NIMHD_LR), and the red curve is a zoom-out, branched out of run NIMHD_LR and run at ℓmax = 23 (run NIMHD_LR_2). The significant overlap between the results of NIMHD_LR and NIMHD_LR_2 (green and red curves) shows that the results of this doubly zoomed-out run are converged with regard to ρacc.

In the text

	Fig. D1. Real part of the Toomre Q parameter of the circumstellar disk of run NIMHD, seen in a top-down view where each panel corresponds to a different time, with t = 0 corresponding to the epoch of protostellar birth. Only cells belonging to the circumstellar disk were selected when producing these plots.
In the text

	Fig. E1. Average radial mass flux within the circumstellar disk in run NIMHD, measured in radial bins at t ≈ 0.5 yr (where t = 0 corresponds to the epoch of protostellar birth). Measurements are performed along the upper layers of the disk (solid line) and the main body (dashed line). Only cells belonging to the circumstellar disk were used in the computation.
In the text