© ESO 2021

1 Introduction

It is commonly accepted that the planets of our Solar System (S.S.) formed within 1–100 million yr (Myr), but several lines of evidence show that planetary formation processes may have started on a much shorter timescale. The calcium-aluminum-rich inclusions (CAIs) that define the time origin of the S.S. most likely have formed in a hot environment close to the Sun. To explain their presence in carbonaceous chondrites far from the Sun, it has been proposed that they formed in the disk during the collapse stage of the prestellar core, when infall of diffuse medium was significant. They were subsequently transported outward during the viscous expansion of the disk (see e.g. Yang & Ciesla 2012; Desch et al. 2018) or were transported in disk winds or stellar outflow (Shu et al. 2004; Yang & Ciesla 2012; Pignatale et al. 2018). Moreover, simulations of stellar cluster formation within a molecular cloud showed that the typical infall time onto individual star-disk systems is about a fraction of a Myr (Padoan et al. 2014).

The low abundance of water in the inner S.S. (Morbidelli et al. 2016; Hyodo et al. 2019) and the identification of several isotopic reservoirs in the S.S. (Kruijer & Kleine 2017; Van Kooten et al. 2016) argue for an early formation of a dynamical barrier such as Jupiter at about one Myr, but the building blocks of a barrier like this must have formed much earlier. Finally, a recent study showed that disks do not contain enough detectable dust to form the exoplanet populations (Manara et al. 2018), implying that either (1) most planets form in a few 100 thousand years (kyr) at the class 0 or class I phase, or (2) a continuous dust replenishment process feeds planets over the lifetime of the disk (Manara et al. 2018). All these elements raise the fundamental question whether planetary accretion started during the assemblage of the protoplanetary disks or after, as is often assumed. If accretion did start in the disk during the prestellar core infall, this would be a major change in our understanding of planet formation, as most studies still presume initial conditions from the Minimum Mass Solar Nebula model (MMSN, Hayashi et al. 1979; Weidenschilling 1977), which is pertinent only for an isolated disk.

2 Current understandings of the protoplanetary disk and the scope of the present work

2.1 Beyond the alpha-disk model

Protoplanetary disks are classically described using the alpha-disk model (Shakura & Sunyaev 1973), which considers an isolated disk. While this model is highly idealized, it is a useful tool for studying planet formation because it is flexible, and it is still the basis of most dustgrowth studies (see, e.g., Birnstiel et al. 2012; Testi et al. 2014, for a review) and planet population synthesis relevant for exoplanets (Alibert et al. 2013). Accretion and transport processes across the disk are described with the effective viscosity $\nu =\alpha c_{\text{s}}^2/\Omega$, where α is the dimensionless local shear stress coefficient that we define below, c_s is the thermal sound speed, and Ω is the Keplerian frequency.

The accretion rate onto the star in steady state, Ṁ = 3πνΣ, is also controlled by α. The disk is heated by both the protostellar radiation and the viscous dissipation (Hueso & Guillot 2005). The value of α is linked to the intensity and the structure of turbulence in the disk. It has been thought during the last 15 yr that the source of turbulence is magnetorotational instability (MRI), giving α > 10⁻³ (see, e.g., Balbus & Hawley 1991; Fromang & Nelson 2006) in ideal magneto-hydrodynamic (MHD) cases. However, later it was realized that nonideal MHD effects (ambipolar diffusion, ohmic diffusion) may preclude disk magnetization, such that α may in fact be ≪ 10⁻³ in the planet-forming region (see, e.g., Turner et al. 2014, for a review). A low turbulence state like this may favor dust coagulation and planetesimal formation (Bai & Stone 2010; Youdin 2011), but a low value of α cannot account for observed accretion rates onto the star, especially for young disks, with 10⁻⁸ M_⊙ yr⁻¹ < Ṁ < 10⁻⁵M_⊙ yr⁻¹ (Hartmann et al. 1998).

Studies of disks with nonideal MHD show that the disk may in fact be vertically stratified, with a low-α midplane, high-α upper layers that drive accretion (Turner & Sano 2008; Turner et al. 2014), and possibly additional MHD structures above a few pressure scale-heights (H) such as disk winds (see, e.g., Bai & Stone 2013; Riols et al. 2016). These new understandings lead to descriptions that are more complex than the α-prescription and might be time and position dependent (Suzuki et al. 2016). Nevertheless, all these studies still consider isolated disks.

2.2 Attempts to study an assembling protoplanetary disk

While many works have addressed the issue of disk formation, only a handful have attempted to study early planetary formation processes (dust coagulation, transport, and planetesimal formation) during the collapse of the prestellar core that feeds the growing disk. The main difficulties are (1) modeling the infall of the envelope and (2) connecting this properly to the circumstellar disk. Classically, the self-similar singular isothermal sphere collapse solution has been used (SIS, Shu 1977; Ulrich 1976). Angular momentum conservation inside the system is always assumed (a critical approximation that neglects magnetic forces). Several works showed (Nakamoto & Nakagawa 1994; Hueso & Guillot 2005) that under this hypothesis, disks may form with characteristics similar to the MMSN (Nakamoto & Nakagawa 1994) or to some disks observed in young stellar clusters (Hueso & Guillot 2005). Using the SIS infall solution, it was shown that CAIs could form close to the star from condensing gas and would then be transported away through viscous disk relaxation (Yang & Ciesla 2012; Pignatale et al. 2018).

Drążkowska & Dullemond (2018) showed that planetesimals may form during the gas infall. In these models, the infall from the envelope is represented as a source term for the disk surface density, and the disk transport processes (controlled by α) are set independently. For example, Hueso & Guillot (2005) considered a fully active disk, whereas Pignatale et al. (2018) considered a disk with a dead zone. However, these two critical aspects are in fact linked. For example, as the gas falls onto the disk, it changes the ionization state of the disk and thus α. Second, the turbulent motions inside the disk may be inherited from the turbulence or from the asymmetry in the infalling envelope. In addition, MHD effects may cause some outflow, and part of the gas may be driven above the midplane. Finally, the accretion shock of the infalling envelope may heat the upper layers of the disk, thus changing the thermal structure of the disk. Therefore it is indispensable to properly study the gas infall that is coupled to the transport properties of the disk and the resulting disk structure in terms of the surface density, temperature profile, and α.

2.3 Previous numerical studies of the disk properties that include the envelope effect

Vorobyov & Basu (2010, 2015) studied 2D thin-disk models with radiative cooling, mechanical heating, and irradiation from the protostar and the background. For assumed frozen-in fields, the magnetic tensor was treated as a dilution factor of gravity, and the magnetic pressure was treated as being a few times the thermal pressure. They showed that the accretion processes within the disk are intrinsically variable, while external infall and radiation affect the fragmentation within the disk. Vorobyov et al. (2015) studied the effect of the environment on an embedded disk. The size of the disk is highly subject to the properties of the infalling material, in particular the angular momentum. The disk fragmentation is affected as a consequence. The global mass accretion history onto the star is consistent with that onto the disk, while the burst modes are largely regulated by fragmentation within the disk. This 2D model does not describe the vertical stratification and thus does not have a source term for the surface density. Whether mass transits through the entire disk before reaching the star or falls through upper layers cannot be distinguished in such setup.

Kuffmeier et al. (2017) have developed a zoom-in technique to follow the evolution of nine sink particles at 2 AU in a simulation of a giant molecular cloud with a resolution of 126 AU with an ideal MHD prescription. They showed that the sink mass accretion is highly sensitive to the large-scale environment and is not as simple as an SIS collapse. Kuffmeier et al. (2018) further increased the resolution of one of the zooms to 0.06 AU and demonstrated that numerical convergence is probably not reached. The short-term variation (bursts) of the mass accretion is highly attenuated with increased resolution.

These works have demonstrated the challenges in studying protoplanetary disk evolution during the embedded phase. Compromises including thin-disk approximation and ideal MHD were made to alleviate computational cost.

2.4 Purpose and structure of the present work

The purpose of the present paper is to revisit the classical pure hydrodynamical model of disk assemblage (Nakamoto & Nakagawa 1994; Hueso & Guillot 2005) with 3D nonideal MHD simulations during the first 100 kyr of the building phase of the disk. We wish to investigate (1) the mass flux onto the disk as a function of position and time, (2) the mass accretion rate from the core onto the disk and that from the disk onto the star, (3) whether the turbulent stress inside the disk or the envelope infall rate controls the accretion rate onto the star, (4) the lost of angular momentum during the infall and the accretion, (5) the disk size and its surface density, (6) the temperature profile within the disk and the location of the snow line, and (7) the behavior of α inside the disk.

To avoid confusion, we define “infall” as mass falling from the envelope to the disk, and “accretion” as mass falling from the disk onto the star throughout this manuscript. The difficulty in this type of study is that long-term evolution and detailed structures need to be followed simultaneously, which is a computational challenge. The present work is an attempt in this direction: We present a detailed study of the disk and inflow structure of an assembling protoplanetary disk. The paper is structured as follows: the numerical simulations are presented in Sect. 3. In Sect. 4 we briefly review the classical disk models that have been applied in planetary science. Detailed analyses of the simulation results are presented in Sect. 5. The results and comparison with the α-disk model are discussed in Sect. 6. Comparison to other numerical works is discussed in Sect. 7. We conclude the paper in Sect. 8.

3 Numerical simulations

To study the formation of the protoplanetary disk self-consistently, we simulated the collapse on the scale of a prestellar core. The code RAMSES (Teyssier 2002; Fromang et al. 2006), which treats MHD with an adaptive mesh refinement (AMR) scheme, was used. Because the initial prestellar core is of subparsec size and the disk forms as a high-density region of a few dozen AU, this collapse problem implies a wide dynamical range of scales, equivalent to four orders of magnitude. This type of simulation has been proven to be computationally difficult because the required resolution is high, which means both high computational power and large memory. Even with an AMR scheme, the highly structured mesh configuration makes the parallelization of the computation inefficient. For this reason, the simulations in this study were typically performed on about 100 cores at most.

3.1 Physical setup

The initial condition was a prestellar core with a flat density profile, 0.02 pc (3800 AU) radius, and total mass of 1 M_⊙. The core initially had a weak solid-body rotation around its center. The radial profile of the rotation is not a major concern because the initial amount of rotation does not strongly affect the outcome of the disk (Hennebelle et al. 2020a, Paper I hereafter). A uniform magnetic field threaded the core and was inclined by θ_mag = 10° with respect to the rotation axis. The magnetic field strength is described with the mass-to-flux ratio, μ, which is the normalized value with respect to its critical value. The value μ = 1 implies magnetic criticality, and an increase in mass or decrease in field flux leads to gravitational collapse. The choice of 1∕μ = 0.3, compatible with observations (e.g., Crutcher et al. 2004; Maury et al. 2018), gives a field of 0.2 mG at the center of the core. This run had no initial turbulent velocity field, although turbulent motion was generated spontaneously during the collapse as a result of density fluctuations.

The physical properties were defined with dimensionless numbers. The thermal virial parameter α_vir = 0.4 specifies the ratio of thermal to gravitational energy, and β_rot = 0.04 is the ratio of rotational to gravitational energy. A weak m = 2 density perturbation of 10% was introduced perpendicular to the rotation axis to break the symmetry during the collapse. This prevents artificial symmetric pattern formation due to the grids. The effects of these physical parameters are discussed in Paper I.

We solved the complete MHD equations considering the ambipolar diffusion (ion-neutral friction, Masson et al. 2012, 2016; Marchand et al. 2016). The temperature was calculated with full radiative transfer with a flux-limited diffusion scheme (FLD, Commerçon et al. 2011), and the temperature floor was set at 20 K. The thermal irradiation from the protostar was also considered using the prescription of Kuiper & Yorke (2013). The accretion luminosity was not included here because its actual value is highly uncertain (see the discussion in Hennebelle et al. 2020b).

We did not take ohmic dissipation and Hall effect into account when treating nonideal MHD. Ohmic dissipation is important at higher densities and induces the formation of a small secondary disk around the protostar (Tomida et al. 2015). Including this effect might not affect our general conclusions because the central high-density part is highly demagnetized, as we show in Sect. 5.4.2. On the other hand, the Hall effect seems to indeed play a role during the formation of the disk (Tsukamoto et al. 2015; Marchand et al. 2018; Zhao et al. 2020), and this remains to be considered in future studies. It is hampered by large uncertainties, however, because knowledge on the grain properties is limited.

3.2 Numerical setup

The AMR scheme requires that the local Jeans length is resolved by 20 cells. The base grid had 2⁶ = 64 cells in each dimension, and the canonical run had a maximum refinement level ℓ_max = 14, equivalent to a resolution dx = 0.93 AU. With the canonical setup, the interior of the disk is not very well resolved, and the dynamics within the disk is not fully physical. Therefore the analyses of this run focused more on the infall of the prestellar core envelope onto the disk and on the formation of the disk. This canonical choice for the resolution allowed us to follow the disk evolution during almost 100 kyr (150 kyr after the simulation started). Increasing the resolution requires a much longer computational time, and we therefore only reran the simulation from some snapshots with increased resolution to examine the consequences and to resolve the structure of the disk.

At the center of the collapse, where the high density is no longer resolved by the fluid description, we inserted a sink particle to trap the collapsed mass (Bleuler & Teyssier 2014). This was necessary because otherwise, the warm and dense gas accumulated at the center expands in a way that is not physical. The sink particle forms when the gas molecular number density exceeds n_thres = 3 × 10¹³ cm⁻³. It accretes from its surrounding within 4dx. A thresholdaccretion scheme was used such that every cell having density higher than n_acc = n_thres∕3 gives a fraction of the excessive mass to the sink particle. This fraction, c_acc, was set to 0.1 in this study.

3.3 Restart simulations with increased resolution

In the canonical run, the entire disk is refined at the highest resolution (2¹⁴). The vertical structure of the disk is poorly resolved, however. To study the detailed dynamics of the disk interior, two restarts were performed at about 40 (R_40ky_ℓ18) and 80 (R_80ky_ℓ18) kyr after sink particle formation (at simulation time 61 kyr) with increased resolution. The model parameters are summarized in Table 1.

Figure 1 shows four column density snapshots of the disk with edge-on and face-on views. The first snapshot corresponds to the time immediately after the sink particle formation. Shortly before the sink particle forms, the dense core grows and flattens as a result of the infall of the envelope. At the same time, the accreted gas with angular momentum leads to the formation of a rotating protoplanetary disk around the central dense core. This core-disk structure is initially thick and reaches a mass of almost 0.1 M_⊙. After the central star forms, most of the mass of this dense structure is quickly accreted, and the disk mass drops and stays at around 2− 3% of the stellar mass during the following evolution.

Despite the inclined initial magnetic field, the disk is almost aligned with the prescribed rotation. In the long-term evolution, the disk gradually grows in size from a few dozen to almost 100 AU at simulation time 120 kyr, while there are fluctuations due to the turbulent nature of the infalling envelope. The disk sometimes decreases in size as a result of sudden strong accretions that are often asymmetric, that is, that come from one side of the disk, and as a result of the interchange instability. Magnetic field is also expelled as an outflow in the form of loops or bubbles (see, e.g., the second row of Fig. 1), which is highly dynamic.

Figure 2 shows the restart at 40 kyr along with the canonical run. With increased resolution, the disk structure is better resolved, while the global profile stays similar. The refinement is increased in steps, passing progressively from level 15 to level 18, in order to avoid numerical artifacts from abruptly increased resolution. The final resolution at level 18 equals to 0.06 AU. When we increased the refinement levels, we had to adjust n_acc correspondingly to avoid artifacts. As explained in Paper I, the disk mass is very sensitive to n_acc. Estimated with a simple analytical model, n_acc should scale roughly as dx⁻².

We caution here that the value of n_acc likely sets an inner boundary condition for the disk and determines the density profile throughout the disk, and the temperature is affected in turn (see discussions in Hennebelle et al. 2020a). Machida et al. (2014) and Vorobyov et al. (2019) have also shown similar results. When disk simulations are interpreted physically, we therefore need to be cautious. As we show below, the effect of the choice of n_acc is seen in the disk global density, in the evolution of the snow line, and in the mass accretion rate onto the central star.

4 Existing models of disk formation

The collapseof the prestellar core is usually treated with axisymmetric and unmagnetized assumptions because they are simple and we lack observational evidence to support the contrary. Before analyzing our simulations, we briefly review what has been learned from existing models.

4.1 Collapse of a prestellar core: a purely hydrodynamical scenario

One classical example of the collapse of a prestellar core is the SIS (Shu 1977), which describes an inside-out collapse of mass shells that accrete onto a central star with a time-independent mass accretion rate, $\dot{M}=0.975c_{\text{s}}^3/G$, where G the gravitational constant. While it should be kept in mind that this solution is just one special case among a whole family of self-similar solutions (e.g., Whitworth & Summers 1985), it has been applied to many models for its simplicity. Furthermore, in reality, a prestellar core has a finite extent and the mass accretion rate onto the central object decreases during the collapse and eventually stops, which results in an isolated disk. Observations have measured mass accretion rates a few times higher than this canonical value (~ 2 × 10⁻⁶ M_⊙ yr⁻¹ for isothermal gas at 10 K) for class 0 objects that are at the beginning of their protostellar collapse (e.g., Ohashi et al. 1999; Hirano et al. 2002; Ward-Thompson et al. 2007). Moreover, a declining mass accretion curve has been observed for young stellar objects between class 0 and class I stages, and simple accretion models have been proposed (e.g., Henriksen et al. 1997; Whitworth & Ward-Thompson 2001).

An easy although not unique (see, e.g., Verliat et al. 2020) way to form a disk is to introduce a small amount of rotation into the prestellar core under the assumption that the collapse solution remains unaffected at large scales. Because angular momentum is conserved during the collapse (with at least about 100 times the contrast in size), rotation is significantly amplified, and this leads to the formation of a flattened structure in the centermost region around the star. While some numerical works exist (e.g., Li et al. 2014), very few studies have examined analytically how a prestellar core collapses to form a star that is surrounded by a protoplanetary disk because the breakup of spherical symmetry complicates the problem.

This simplified picture consists of prescribing each collapsing mass shell with a certain angular velocity and deriving how this mass is distributed radially when it arrives in the equatorial plane. The accretion onto the disk is expressed as the source term S(r, t), which describes the mass flux as a function of the radial position and time. We discuss two widely accepted propositions. These models provide in a deterministic manner (1) the total mass accretion rate onto the disk as a function of time, (2) the radius of the accreting zone of the disk as a function of time, and (3) the spatial distribution of mass flux onto the disk surface.

Fig. 1 Column density snapshots of the canonical simulation R_ℓ14. Left: edge-on. Right: face-on. The coordinates are in AU, and the origin is at the box center. The images are centered on the sink particle, and the white dashed circle indicates the numerical radius of sink accretion (4dx). After the formation of the sink particle, the originally roundish structure of the flattened first Larson-core is quickly accreted and a flat disk forms, with sightly smaller size. The protoplanetary disk does not stay at the center of the computation box and might experience a complex accretion history. The size of the disk grows slowly from its initial radius of about 20 AU and presents fluctuations. At about 120 kyr, the disk grows more rapidly and reaches almost 100 AU. Although only a few snapshots are presented here, the actual picture is highly variant in time.

Fig. 2 Zoomed views of the disk at age (103.1 kyr in simulation time) from edge-on (left column) and face-on (right column). Top row: R_ℓ14. Bottom row:R_40ky_ℓ18. The column density is shown as in Fig. 1. The disk interior is clearly more resolved in R_40ky_ℓ18 with 0.06 AU resolution, particularly in the vertical direction.

4.1.1 Free fall along a parabolic trajectory

As proposed by Ulrich (1976), in spherical coordinates (radius r, latitude θ, and longitude ϕ), all particles passing through radial position r₀ are assumed to have angular velocity $\stackrel{.}{{\varphi }_0}$ and move in the plane that is instantaneously defined with the vectors r and ϕ. A parabolic trajectory, with the central star as the focus, is uniquely defined with these two criteria, while the nonzero velocity in the two other directions are not explicitly specified. This spherical shell is not a rotating solid body, but this did not concern the author, probably because the introduced amount of rotation is small. The particle is considered as accreted onto the disk when its trajectory intercepts the equatorial plane. The radial extent, or the centrifugal radius, of this region of accretion is derived as $$r_{\text{c}}=\frac{{\stackrel{.}{{\varphi }_0}}^2{r}_0^4}{GM},$$(1)

which is the pericenter of the particle that arrives from the equatorial plane, with M being the mass of the star-disk system, which is assumed to be dominated by the stellar mass.

When we assume initial uniform distribution of mass on the shell, the source function can be derived as $$S_{\text{p}}(r,t)=\frac{\dot{M}(t)}{4\pi r_{\text{c}}^2(t)}{\left(\frac{r}{r_{\text{c}}(t)}\right)}^{-1}{\left(1-\frac{r}{r_{\text{c}}(t)}\right)}^{-1/2},$$(2)

where Ṁ(t) is the mass accretion rate of the mass that falls from the radius r₀, and is assumed to be a constant of time in most of the existing models, as previously explained.

4.1.2 Free fall onto a Keplerian radius

On the other hand, Nakamoto & Nakagawa (1994) and Hueso & Guillot (2005) assumed that the particle, conserving its angular momentum, arrives at the radius of Keplerian rotation, regardless of the trajectory. With the same setup of initial mass and angular momentum distribution as in the previous model, the source function becomes $$S_{\text{k}}(r,t)=\frac{\dot{M}(t)}{8\pi r_{\text{c}}^2(t)}{\left(\frac{r}{r_{\text{c}}(t)}\right)}^{-3/2}{\left(1-\sqrt{\frac{r}{r_{\text{c}}(t)}}\right)}^{-1/2}.$$(3)

4.1.3 Implication of these formalisms

In the model by Ulrich (1976), the particle has sub-Keplerian rotation when it reaches the disk and will fall inward. The exchange in angular momentum caused by friction of particles in the disk is not described, and a disk evolution model is further required to follow the redistribution of matter inside the disk. In the model by Nakamoto & Nakagawa (1994) and Hueso & Guillot (2005), on the other hand, the trajectory is not described and the particle directly reaches its centrifugal radius, giving a slightly more centrally concentrated source function.

The two models are shown in Fig. 3. For a more intuitive understanding, we also show the accumulated mass accretion within radius r, $$\dot{M}(<r)={\displaystyle {\int }_0^{r}S}(r^{\prime })2\pi r^{\prime }\text{d}{r}^{\prime },$$(4)

which equals zero at the origin and the total Ṁ at r_c. The difference between the two prescriptions does not appear to strongly affect the evolution of the disk (Hueso & Guillot 2005). Both models rely on pure hydrodynamical arguments, and the initial cloud is mapped to the disk deterministically, in a shell-by-shell manner. The instantaneous source function entirely depends on the assumption of uniform mass distribution on the shell and constant angular velocity throughout the shell. While widely accepted because they are convenient, these models are probably too simplistic.

The mass accretion rate Ṁ(t) and centrifugal radius r_c(t) are both determined by the initial conditions (density and angular momentum profiles) of the core. The inside-out collapse solution of Shu (1977) suggests a constant accretion rate Ṁ (t). By assuming a rotation profile of the prestellar core (often solid body, which leads to r_c (t) ∝ t³), the full accretion history of the disk can be constructed. Simple variations of this class of models can be prescribed with different initial density and rotation profiles, leading to different histories of Ṁ (t) and r_c (t).

The source function is of particular interest for the studies of early S.S. evolution because it implies the temperature-pressure history that the matter inside the S.S. can experience. In a purely hydrodynamic model, a non-negligible fraction of mass arrives at the inner part of the disk, where the temperature is high due to viscous heating and protostellar irradiation. The centrifugal radius grows rapidly in time such that accretion later arrives at the cooler part of the disk. This model has been used to explain the composition diversity of chondritic meteorites (Pignatale et al. 2018).

Fig. 3 Upperpanel: source function as proposed by Ulrich (1976) and Hueso & Guillot (2005). The radius is normalized to the centrifugal radius rc, and the source term is normalized to the mean surface mass flux (${\dot{M}}_{\text{tot}}/(\pi r_{\text{c}}^2)$) such that its surface integral gives unity. Lower panel: cumulated mass accretion rate (normalized by the total mass accretion rate), Ṁ, integrated within radius r. The values of Ṁtot and rc are to be specified as functions of time according to the initial conditions of the prestellar core envelope.

4.2 Effects of magnetic field and nonideal MHD

Observational evidence increasingly shows that the star-forming gas is likely magnetized, and so is the disk because it is unlikely that the magnetic field at this scale is completely removed. In the past decade, numerical simulations started to take nonideal MHD effects in the studies of early disk formation into account. The authors found that the size of the disk is regulated by the magnetic field and is smaller than what would have been expected when pure hydrodynamics were considered (see e.g., Dapp & Basu 2010; Inutsuka 2012; Li et al. 2014; Tsukamoto et al. 2015; Tomida et al. 2015; Hennebelle et al. 2016; Machida et al. 2016; Zhao et al. 2018; Wurster & Li 2018; Lam et al. 2019). This is also in line with the observations suggesting that most class 0 disks are small (<50 AU, Andrews et al. 2018; Maury et al. 2019).

In this work, we aim to provide a revised model for the assembly of a protoplanetary disk by taking into account the nonideal MHD effects and the finite extent of the prestellar core. This can serve as more realistic framework for future studies of the early evolution of protoplanetary disks.

5 Disk characteristics evaluated from simulations

Before performing the canonical run, we simulated the collapse of a nonmagnetized core. The outcome is very similar to what is suggested by classical models: a disk forms due to angular momentum conservation, and its size grows quickly in time. However, this disk is highly unstable and fragments very quickly (see, e.g., Matsumoto & Hanawa 2003). Therefore it is unlikely that this scenario canlead to an extended and evolved disk. Observations (e.g., ALMA Partnership 2015) also suggest that disks are probably small at early times. This motivated us to include the magnetic field along with its nonideal effects.

In the presence of ambipolar diffusion, the disk size is regulated by the competition between the magnetic braking that decreases the angular momentum of the gas and the ambipolar diffusion that leaks the field lines outward and reduces the braking (Hennebelle et al. 2016). Given the very different picture with respect to ideal MHD simulations, nonideal MHD effects are indeed necessary to understand the evolution of disk structure. The magnetic flux is allowed to leak outward, reducing the braking, and a rotation-supported demagnetized disk can form as a consequence.

The simulation started with a core in free collapse, and thus it takes some time for the mass to accumulate before the protostar and the disk form. When the central density increases, the gas is heated by the collapse due to the increasing dust opacity, and the pressure locally supports against the self-gravitating collapse. A quasi-stationary first hydrostatic dense core (or first Larson core, Larson 1969; Masunaga et al. 1998; Vaytet & Haugbølle 2017) forms as a consequence, and this core has a slightly flattened shape because of rotation. Further mass accumulation leads to the formation of a sink particle that represents the protostar at the center about 61 kyr after the beginning of the simulation.

We evaluated several disk properties from the simulations. All quantities were averaged in the azimuthal direction and between the upper and lower halves of the disk.

5.1 Disk geometry and global properties

The first step was to clearly identify the disk region and infer its geometrical properties. For this purpose, we analyzed the distribution of mass as a function of radius and altitude to derive the scale height, surface density, and radial extent of the disk. The analyses were made inside a selected cylindrical region with radius r_cyl = 100 AU and half-height z_cyl = 100 AU, of which the center and the axis were first evaluated by selecting dense cells over a certain threshold as a first step of disk identification (see details in Appendix A). At later times when the disk had grown in size, the selected region was increased to r_cyl = z_cyl = 1000 AU. As we show below, because the disk is mostly dominated by its central dense part, the size of the cylindrical region is not important as long as it is large enough with respect to the disk.

5.1.1 Vertical density profiling: scale height H

We calculated the azimuthally averaged vertical density profile for all radii. The binning size in the radial direction was taken to be dx, the smallest cell size (or several times dx if the disk was large with respect to the resolution). We evaluated the disk scale height with several approaches, which were all expected to give identical results if the disk is vertically isothermal and follows an exact Gaussian density profile in pressure equilibrium. The thermal scale height from vertical hydrostatic equilibrium (assuming a dominating stellar mass) is $$h_{\text{p}}(r)\approx \sqrt{\frac{c_{\text{s}}^2(r)r^3}{G[M_{\ast }+M_{\text{d}}(r)]}},$$(5)

where c_s(r) is the local sound speed, r is the cylindrical radius, M_* is the stellar (sink) mass, and M_d(r) is the disk mass inside r. The disk is not a point mass, and the complete gravitational field should be calculated by integrating over the whole disk. We used this approximate formula for simplicity because the disk mass is only a few percent of the stellar mass. Because the disk is not exactly isothermal, the thermal sound speed was calculated as $c_{\text{s}}^2={\displaystyle \sum P\text{d}V}/{\displaystyle \sum \text{d}m}$, where P, dV, dm are the pressure, volume, and mass of each cell, and the summation was performed vertically within cylindrical shells. Here again, we verified thatthe vertical upper limit of the sum has no significant effect, and we simply summed over the whole cylinder.

In the case of a vertical Gaussian profile, $$\rho (z)=\frac{{\Sigma }_{\text{tot}}}{\sqrt{2\pi }H}\exp \left[\frac{-z^2}{2H^2}\right]$$(6)

the scale height can be recovered from the first or second moment of mass: $$H=\sqrt{\frac{\pi }{2}}\frac{1}{{\Sigma }_{\text{tot}}}{\displaystyle {\int }_{-\infty }^{\infty }\rho }z\text{d}z={\left[\frac{1}{{\Sigma }_{\text{tot}}}{\displaystyle {\int }_{-\infty }^{\infty }\rho }{z}^2\text{d}z\right]}^{1/2}.$$(7)

In practice, we calculated $$h_1(r)=\sqrt{\frac{\pi }{2}}\frac{{\displaystyle \sum z\text{d}m}}{{\displaystyle \sum \text{d}m}} \text{and} h_2(r)={\left[\frac{{\displaystyle \sum z^2\text{d}m}}{{\displaystyle \sum \text{d}m}}\right]}^{1/2},$$(8)

where the summation was performed vertically within cylindrical shells. These quantities were used as estimates of the thermal scale height.

Figure 4 shows snapshots of the scale height calculated for the two high-resolution restarts compared to the canonical run. These results of h_p, h₁, and h₂ indeed reasonably coincide with one another inside the disk, and the significant incoherence marks the radius at which the disk ends (left dotted vertical line). This also confirms that our choice of the cylindrical region size is sufficient to avoid errors because the Gaussian density profile drops very quickly at high altitudes.

The two vertical dotted lines denote the disk characteristic radii that we define in Sect. 5.1.2 and further discuss in Sect. 5.1.3. The first radius, R_kep, corresponds to the inner part of the disk that is centrifugally supported in the radial direction and thermally supported in the vertical direction (see Sect. 5.4.1). The scale heights determined from three different definitions thus coincide very well with one another. The second radius, R_mag, corresponds to an outer region of the disk where the magnetic support is relatively strong, and deviations start to appear. Beyond R_mag, the moments of mass no longer give much information as the infalling envelope does not have a flattened geometry. For a vertically uniform density profile (which is almost the case in the outer part of the disk and the envelope), $h_1=\sqrt{\pi /8}{z}_{\text{cyl}}\approx 0.577z_{\text{cyl}}$ and $h_2=\sqrt{1/3}{z}_{\text{cy}}\approx 0.626z_{\text{cyl}}$. They are almost indistinguishable either inside or outside the disk, thus their discrepancy with h_p exhibits the clear change to a nonthermally supported region.

The overall behavior is similar in the high- and low-resolution runs. However, the high-resolution runs provide more precise descriptions of the disk height profile and allowed us to probe the high-density interior of the disk. The horizontal gray line in Fig. 4 marks the resolution dx. The disk in the canonical run contains only one cell in the vertical direction within a radius of 10 AU, while the high-resolution restarts resolve the disk down to radius of 1 AU. Run R_40ky_ℓ18 shows a larger and better delineated transition zone between the two characteristic radii, while in the canonical run this region is barely resolved by a few cells. The disk has nearly constant aspect ratio H∕R ≈ 0.1 up to almost 20 AU radius. On the other hand, the disk has a larger radius after simulation time 120 kyr, and thus the run R_80ky_ℓ18 resolves the disk center better, but does not give much additional information at large radius because it is already well resolved in the canonical run.

5.1.2 Radial density profiling: surface density Σ

The surface density of the disk is obtained as $${\Sigma }_{\text{cyl}}(r)=\frac{{\displaystyle \sum \text{d}m}}{2\pi r\text{d}r}={\Sigma }_{\text{tot}}(r) \text{erf}\left(\frac{z_{\text{cyl}}}{\sqrt{2}\text{H}}\right)\approx {\Sigma }_{\text{tot}}(r),$$(9)

where Σ_cyl is the column density within the cylindrical region, and Σ_tot is the total surface density if the vertical density profile followed an exact Gaussian, with the mass contained within ± z_cyl equaling the numerically measured value. The error function is a correction for the mass outside the cylinder vertical extent that is not taken into account. However, this correction is almost negligible for H ≪ z_cyl.

Figure 5 shows snapshots of the Σ profile. At 138 kyr, the surface density profile is close to a power law with a sharp drop near 13 AU, and it becomes flat again outside 30 AU. The rapidly decreasing part is due to increased magnetic support toward the outer disk (see Sect. 5.4.2 for discussions of the magnetic field). Most of the disk mass is within the inner region, while sometimes the region with sharply decreasing surface density can contain up to a quarter of the total disk mass. The outer flat part, on the other hand, is not a real surface density, but a truncated value of the diffuse envelope that is more vertically extended. We performed a three-segment power-law fit to the measured Σ_cyl(r), and the transition points define the two characteristic radii, R_kep and R_mag, which we discussed in the following paragraph.

Fig. 4 Scale heights from vertical hydrostatic equilibrium, hp (blue), from the first (h1, orange) and second (h2, green) moment of mass, plotted against the radius. Top two panels: R_ℓ14 and R_40ky_ℓ18 at 103 kyr. Bottom two panels: R_ℓ14 and R_80ky_ℓ18 at 138 kyr. The horizontal and the left vertical gray lines mark the numerical resolution (dx, outside the figure extent, thus not shown for the high-resolution runs). The values are averaged within rings of thickness dx and plotted against the mean radius of the rings. The second vertical gray line marks the radius of the numerical sink accretion scheme (4dx). The two vertical dotted lines mark the characteristic radii of the disk (Rkep and Rmag), as we defined in Sect. 5.1.2. Within Rkep the disk is close to vertical thermal equilibrium and the three definitions of scale height coincide very well. Beyond this radius, magnetic support becomes strong and deviation occurs. The thermal scale height, hp (r), is fitted to apower law, which is shown in the legend (in units of AU).

5.1.3 Disk radius R

Two characteristic radii, R_kep and R_mag, were definedfrom a three-segment power-law fit of the surface density profile (Sect. 5.1.2). This is shown in Fig. 5 with the vertical dotted lines. The inner part of the disk has a shallow power-law profile and is almost nonmagnetized. The outer part, on the other hand, is supported by the magnetic field, and the surface density drops quickly (see Sect. 5.4.2). These two characteristic radii are marked in all the following radial profile figures to show that they indeed correspond to some transitions of the properties inside the disk. We recall that no clear-cut disk radius exists, whileour choice of characteristic values allows making reasonable evaluations of other disk properties in our analyses.

5.2 Global evolution of the disk: mMass and radius

With the nonideal MHD effects, the disk becomes self-regulated and its size evolves slowly. This is due to the balance between the rotation and the magnetic braking as well as to that between the magnetic field induction through rotation and its loss through diffusion (Hennebelle et al. 2016). The radius very quickly reaches ~ 20 AU and becomes quasi-stationary. The disk receives mass from the infalling envelope, and at the same time, some of its mass is lost through accretion onto the central star. The mass inside the protoplanetary disk stays at ~ 0.01 M_⊙.

Figure 6 shows the disk radius evolution for the canonical simulation at 1 AU resolution, with R_kep plotted as a solid blue line and R_mag as a dashed blue line. The sink particle forms at simulation time 61 kyr, and we only show measurements from ~75 kyr, when most mass of the star-disk system is accreted by the sink particle and a flattened disk is well defined. Figure 1 shows that the star-disk system is more like a flattened ellipsoid at 61 kyr when the sink particle has just formed. The size, R_kep in particular, does not vary strongly within more than 50 kyr. At about 120 kyr, the disk starts to grow in size up to ~50 AU. The magnetized part grows significantly, reaching several times the size of the Keplerian disk. The growth of the disk is possibly a consequenceof the evolution of the surrounding medium properties, and in particular, of the magnetic field. It may also indicate that the coupling with the surrounding envelope, with which angular momentum is exchanged, becomes less tight, and thus the efficiency of magnetic braking is reduced. This particular point certainly requires further investigation. The high-resolution restarts are overplotted and the disk sizes are consistent. With increased resolution, the time steps are significantly reduced such that the evolution is not followed during very long time spans.

Figure 7 shows the mass evolution of the star, the disk mass within R_kep, and that within R_mag. Before the sink formation, a flattened dense structure is formed. Quickly after sink particle formation, most of the mass is accreted by the star and the disk mass drops to stably ~ 2−3% of the stellar mass. Most of the mass is contained within R_kep, where the disk is essentially Keplerian. When the disk starts to expand at later times, the magnetized region grows significantly, and its masssometimes reaches one-third of the mass of the Keplerian part. The overplotted high-resolution runs are essentially identical.

Figure 8 shows the mass accretion rate of the sink particle and that across the disk surface (at min (3h_p, 4dx)). The accretion rate of the sink is directly derived from the sink mass history. The occasional bursty behavior is probably due to the mass accumulating at the inner disk boundary before the density exceeds n_acc. Kuffmeier et al. (2018) also saw a similar behavior in their zoomed-in systems at 2 AU resolution, but this feature diminished with increased resolution. Therefore this phenomenon should be interpreted with care. The two values are broadly comparable. The mass accretion onto the disk surface slowly decreases in time. On the other hand, the sink accretion rate experiences a more serious drop. This is probably because the accretion becomes numerically more difficult as the disk density slowly decreases while n_acc is kept a constant. At later times, Ṁ_Disk > Ṁ_sink. This is also reflected in the growth of mass and size of the disk (Figs. 6 and 7).

The mass accretion rate across the disk surface in the high-resolution runs is compatible with the canonical run. However, R_40ky_ℓ18 presents zero sink accretion. This reflects the sensitivity of the sink evolution to the numerical parameters, and n_acc was probably not very well adjusted for this case (see Sect. 3.3). On the other hand, R_80ky_ℓ18 shows a higher sink accretion rate than the canonical run. This might be due to the better description of the central high-density region, or it might suggest that the value of n_acc is slightly too low.

Fig. 5 Disk surface density profile Σcyl(r) measured inside a selected cylindrical region (blue). Top two panels: R_ℓ14 and R_40ky_ℓ18 at 103 kyr. Middle: radially accumulated mass of the disk of R_40ky_ℓ18 at 103 kyr. Bottom two panels: R_ℓ14 and R_80ky_ℓ18 at 138kyr. A three-segment power-law fit is overplotted (orange), and the power-law exponentsare shown in the legend. Two vertical dotted lines show the characteristic radii Rkep and Rmag that correspond to the transition of the power-law slope. The masses of the star M* and of the disk Md (inside Rkep plus betweenRkep and Rmag) are also displayed. Most of the disk mass is contained within Rkep.

Fig. 6 Radius evolution of the disk. The radii measured from simulations as described in Sect. 5.1.3 are plotted with solid (Rkep) and dashed (Rmag) lines for R_ℓ14 in blue. The disk initially has a Keplerian radius of about 20 AU for almost 50 kyr and a magnetized envelope that is roughly 1.5 times larger. At later times, the disk grows to nearly 50 AU in radius, and the magnetized region grows even more significantly to more than twice the size of the demagnetized Keplerian disk. Runs R_40ky_ℓ18 and R_80ky_ℓ18 are overplotted (with the same line styles), and the size of the disk is consistent with different resolutions.

Fig. 7 Mass evolution of the star-disk system. Run R_ℓ14 is plotted as a thick blue line for the sink mass, a thin line for the mass within Rkep, and as a dashed line for the mass within Rmag. A flattened core-disk structure forms at the beginning. After the sink forms, it quickly accretes most of the mass inside the disk, and the disk mass drops to about 1% of the stellar mass. The magnetized region around the Keplerian disk, though significant in size, contains much less mass. The high-resolution restarts are overplotted with the same line styles, and there is an overall good agreement.

Fig. 8 Mass accretion rate evolution. The canonical run at 1 AU resolution is plotted in blue. The mass accretion rate of the sink particle (thick line) is derived directly from its mass evolution. The accretion rate onto the disk is the integral of the measured source function (dashed line, as described in Sect. 5.5.3). The mass accretion rate is initially high an decreases slowly with time. The mass accretion onto the sink particle, though generally decreasing in time, shows some bursty behavior. This is due to the accumulation at the inner disk border shortly before reaching the accretion threshold. After 130 kyr, the density of the disk drops so much such that the threshold density of the sink accretion algorithm is no longer met and the sink hardly accretes. The high-resolution restarts are overplotted (solid line for sink accretion and stars for disk surface accretion). The mass flux onto the disk surface is unchanged in general, while the sink accretion is sensitive to the inner boundary condition of the disk (i.e., the choice of nacc, see Sect. 3.3).

5.3 Temperature inside the disk: snow line and CAI formation

One of the important properties of the protoplanetary disk is the temperature distribution because it is important for the condensation of different minerals, and also for water. Knowing the temperature conditions will allow us to study the formation of precursors of different mineral components of the building blocks of the S.S. The water content, on the other hand, plays crucial roles not only in the conditions of chemical reactions, but also but also for the growth of solid particles because water ice represent the most massive reservoir of solid material in the S.S. In this section, we discuss the thermal structure of the disk and its evolution.

5.3.1 Radial temperature profile

Figure 9 shows snapshots of the radial temperature profile. We show profiles averaged (mass-weighted) within one and three times the disk scale height H, and they are indeed almost identical, confirming the vertical isothermal simplification. There is no properly defined midplane temperature because we only have access to the temperature at the resolution limit, that is, within one cell height. The value of h_p is used for H, thus it is more reliable within R_kep. The gas is relatively diffuse outside and the vertical temperature variation is considerably low. As a verification of the thermal sound speed evaluated in Sect. 5.1.1, we examined the vertical temperature profiles. The temperature only deviates significantly at rather high altitude compared to the disk scale height, and we refer to Appendix B for more details.

In general, the temperature decreases with increasing distance to the central star. This is mostly due to the heating from the central star because the temperature of the disk is significantly lower when the stellar irradiation is not considered (see Figs. 10 and 11). We caution that the kink seen near 4 AU in the canonical run is not physical. It corresponds to the sink accretion radius, where gas behavior is no longer physical. The temperature profile is close to what would be expected for a flared disk illuminated by the star, where T∝ r^−1∕2 (Kenyon & Hartmann 1987). This temperature profile is derived by balancing the flux received at the disk surface (∝ r⁻²) with the local black body radiation (∝ T⁴). More specifically, Chiang & Goldreich (1997) derived $$T\propto r^{-1/2}{\alpha }_{\text{g}}^{1/4},$$(10)

where α_g = dh_p∕dr − h_p∕r is the grazing angle at which the stellar irradiation strikes the disk. Combined with our measurements of h_p ∝ r^1.1−1.3 (see Fig. 4), we can derive a temperature profile slightly shallower than r^−1∕2. However, our disk is surrounded by an envelope that is not optically thin, and the geometry of the disk is not very flared (see Fig. 4). The photons emitted from the protostar are probably unable to penetrate the upper layers freely and directly strike the disk surface at large radii. It is thus more appropriate to describe the temperature distribution in the diffusion regime, where the radial temperature gradient is established through the equilibrium of radiative diffusion, such that (cf. Eq. (5) of Hennebelle et al. 2020b) $$-\frac{4\pi r^2\sigma }{3\kappa (T)\rho (r)}\frac{\partial T^4}{\partial r}=L_{\ast },$$(11)

where L_* is the luminosity from the protostar, κ is the opacity, and σ is the Stefan-Boltzmann constant. The density profile ρ(r) is slightly shallower than r⁻² in general (see Fig. 3 of Hennebelle et al. 2020a). For roughly constant opacity (Semenov et al. 2003) and ρ ~ r^−2~−1, we can derive T ∝ r^{−3∕4~−1∕2}. In conclusion, the disk temperature profile is similar to what would be expected for an irradiated disk, while the physical origin of this profile is different.

For a passive disk that is heated only by viscous dissipation, Bitsch et al. (2013) derived T ∝r^{−3∕4−s∕2} (constant viscosity) or T ∝ r^{−1∕2−2s∕3} (α viscosity), where s is the power-law exponent of the surface density profile Σ ∝ r^−s, depending on different viscosity models. Again, these profiles are very similar to those we derived above. This shows that the temperatureprofile is not a very sensitive proxy with which we can distinguish different heating mechanisms of the disk.

Here weare particularly interested in two temperature values: the water snow line at 150 K, and the CAI condensation temperatureat 1500 K. We used 1500 K as a representative value, while it should be kept in mind that the CAIs are a family of refractory minerals and form in a temperature range around this value. We identified the two radii, R₁₅₀ ~ 7 AU and R₁₅₀₀ < 1 AU. Only in the high-resolution runs (dx = 0.12 AU) is R₁₅₀₀ marginally resolved.

Fig. 9 Radial temperature profiles of the disk. The same snapshots as in Fig. 4 are shown. The temperature is averaged (mass-weighted) vertically within H (orange) and3H (green). The snow-line (R150) and the CAI condensation line (R1500, if present) are marked with vertical dashed cyan lines. In the low-resolution run, the snow-line R150 is very close to the resolution limit and thus not well resolved. High-resolution restarts are indeed necessary to study the interior structure of the disk correctly. The CAI condensation line R1500 appears in R_40ky_ℓ18, while in R_80ky_ℓ18, this line hasmigrated inward due to the decreased overall density and temperature.

Fig. 10 Radial temperature profiles of the disk with and without stellar irradiation at time 120 kyr. Values of R_ℓ14_nofeed are plotted with lighter colors. The temperature is averaged (mass-weighted) vertically within H (orange) and3H (green). The snow line (R150) is marked with a vertical dashed cyan line. The disk temperature is significantly lower when no stellar irradiation is introduced, while the slope of radial profile does not seem to differ significantly.

Fig. 11 Evolution of the water snow line, R150. The snow line is initially located at about 10 AU and migrates inward in time as the disk evolves and loses its mass, finally stabilizing at ~5 AU. High-resolution runs are overplotted. R_80ky_ℓ18 has consistent temperature, while in R_40ky_ℓ18 the snow line migrates to a larger radius, possibly due to the mass accumulation in the disk and the increase in opacity. This reflects the importance of properly choosing nacc with increased resolution (see Sect. 3.3). The run with no stellar irradiation is also plotted for comparison, where the temperature is significantly lower.

5.3.2 Evolution of snow lines

Figure 11 shows the evolution of R₁₅₀, the water snow line. We do not discuss R₁₅₀₀ because it isonly marginally resolved in the runs with the highest resolution.

R₁₅₀ starts at about 10 AU and decreases slowly in time, as the global density of the disk decreases due to accretion onto the star and redistribution of material in the envelope surrounding the disk. After about 80 kyr, the snow line stabilizes at about 5 AU, close to the current location of Jupiter, which may have consequences for its formation (see Sect. 6). R₁₅₀ from high-resolution restarts are overplotted. Run R_80ky_ℓ18 shows a consistent snow-line position. On the other hand, the R_40ky_ℓ18 is significantly heated, and the snow line moves very quickly outward due to the high opacity resulting from the accumulated mass. This reflects that the choice of n_acc is not entirely consistent when we increase the resolution (see Sect. 3.3). We recall that the sink does not accrete in this high-resolution restart (Sect. 5.2). The slight increase in surface density, which is barely perceptible in the Σ profiles, leads to an increased opacity and thus to an increased temperature. Due to the same issue, the evolution trend of the snow line we obtained should be regarded as tentative, and its exact position depends on the detailed physics of the stellar accretion.

To evaluate the significance of the stellar irradiation, we performed run R_ℓ14_nofeed, where the intrinsic luminosity from the protostar was turned off at simulation time 119 kyr. As shown in Fig. 11, the temperature of the disk immediately is significantly lower. This implies that in our simulations, the disk temperature is primarily set by the stellar irradiation. However, because this run was only performed for a short time span, it is unclear whether the disk will reach a different equilibrium state with slightly higher temperature by increasing the density. The effect of irradiation should be investigated in detail, along with the luminosity resulting from mass accretion. We leave this for future studies.

Fig. 12 Velocity profiles of the disk at two resolutions: rotational velocity within H (orange) /3H (green) and infall velocity within H (red)/3H (violet). Both positive and negative values with |u| > 102 are plotted in logarithmic scale, and this is connected with a linear range in between. The Keplerian rotational velocity is plotted in blue for reference. The inner disk is very close to Keplerian rotation, and the deviation becomes remarkable in the outer magnetized part of the disk. The envelope is far from Keplerian rotation, while the infalling gas has roughly half the free-fall velocity. Below the disk accretion radius (vertical gray line), the velocity is not physically meaningful. The disk characteristic radii Rkep and Rmag (the dotted vertical lines) mark the transition between different behaviors.

5.4 Disk dynamical properties

Some disk properties, such as the growth and radial drift of dust and pebbles, as well as planetesimal formation, are closely linked to its dynamics and evolution. Because it is always assumed that the protoplanetary disk has close-to-Keplerian rotation, we first confirm the validity of this common assumption. We then continue to discuss the magnetization and the transport parameters that are essential in controlling dust-growth and transport within the disk.

The azimuthally averaged quantities are presented, and we show both the average within one and three times the disk scale height H. Because 68.3 and 99.7% of the mass are contained within these extents for a vertical Gaussian profile, these quantities should be representative of the disk. All quantities were calculated throughout the selected radius, r_cyl, while the notion of scale height is only valid within the disk radius R_kep. The quantities inferred in the envelope region are therefore only suggestive of the values around midplane.

5.4.1 Velocity: rotation and accretion

We show therotation velocity, u_ϕ, and the radial velocity (parallel to the disk plane), u_r, in Fig. 12. The Keplerian rotation velocity is plotted for reference: $$u_{\text{kep}}(r)\approx \sqrt{\frac{G[M_{\ast }+M_{\text{d}}(r)]}{r}}.$$(12)

This approximation is valid for a disk mass dominated by the stellar mass.

The mass-weighted average is calculated within one and three disk scale heights. The azimuthal velocity is plotted in orange (H) and green (3H). Different behaviors are clearly seen in regions separated by the disk characteristic radii R_kep and R_mag (see the definition in Sect. 5.1.2, and here we see the reason for giving these names). The velocity is very close to the Keplerian rotation between ~ 4 AU (sink accretion radius) and ~13 AU (R_kep). Inside the sink accretion radius, the flow property is not physical and should not be considered. Beyond R_kep, the rotation starts to deviate from the Keplerian value because significant support comes from the magnetic field. Beyond R_mag it is clear that the gas inside the infalling envelope is far from Keplerian rotation.

On the other hand, the infall velocity (negative radial velocity) is plotted in red (H) and purple (3H). In the envelope, the gas falls radially onto the disk at almost half the free-fall velocity, which has the same expression as the Keplerian velocity. In the interior of the disk, the radial velocity is significantly lower than the rotational velocity (by at least two orders of magnitude). The radial velocity flips between inward and outward directions with no clear pattern, and it also varies with time. The difference between the two runs at different resolutions indicates that the internal dynamics is highly variant and has not yet reached numerical convergence. Nonetheless, the innermost disk (below 1 AU) always shows accretion toward the central star. As we discuss in Sect. 5.5.1, most of the mass accretion onto the star occurs directly through a high-altitude channel and does not pass through the bulk of the disk. Therefore the disk is not necessarily always in accretion.

5.4.2 The β parameter: magnetization

The ratio of the thermal and magnetic pressures, $\beta =P_{\text{therm}}/P_{\text{mag}}=8\pi \rho c_{\text{s}}^2/B^2$, indicates which supporting agent is dominant inside and outside the disk region. This is calculated as $$\beta (r)=8\pi \frac{{\displaystyle \sum P}\text{d}V}{{\displaystyle \sum B^2}\text{d}V}.$$(13)

Figure 13 shows β values within H and 3H. Because the temperature inside the disk is high, the thermal pressure largely dominates, with β decreasing outward from 10³ at 1 AU. The value of β crosses unity inside the magnetized zone of the disk (between R_kep and R_mag), marking a transition to the magnetized infalling envelope. As the global collapse proceeds, the central high-densitypart loses the field lines through ambipolar diffusion, and the magnetic pressure is therefore weak in the inner part of the disk. The disk is more magnetized at higher altitudes, giving lower β values at 3H than at H (because the temperature is almost identical to that in Fig. 9). Although significantly demagnetized, the magnetic field is still stronger than what is typically assumed in studies of disk instabilities. For example, Béthune & Latter (2020) used β ∈ [10³, 10⁸].

5.4.3 The α parameter: accretion and diffusion

The α parameter describes the transport in the disk as a result of torques from the turbulence or magnetic field. The cross correlation of different components of a vector field generates a tensor term that implies the transportation of mass and angular momentum (Balbus & Papaloizou 1999). It has contributions from the turbulent fluctuation (Reynolds tensor), the magnetic field (Maxwell tensor), and the self-gravity of the disk. The turbulent diffusion is not very well known within the disk, and thus $\nu ={\alpha }_{\text{Rey}}{c}_{\text{s}}^2/\Omega$ is often used as the turbulent viscosity to describe the diffusion within the disk. We evaluated radial dependence of these terms inside the disk and discuss their effect on the transport of mass and angular momentum.

The stress tensor, $T_{p\varphi }(p,\varphi )\equiv {\langle \rho \delta u_p\delta u_{\varphi }\rangle }_{\varphi }$, is an azimuthal average along a circle perpendicular to the axis, where the subscript p stands for the poloidal direction that has two components r and z, and δ denotes the deviation from the mean value. The dimensionless α expression is defined as $T_{p\varphi }/{\langle \rho c_{\text{s}}^2\rangle }_{\varphi }$, where the denominator is the azimuthally averaged thermal pressure.

The radial mass accretion rate of a stationary disk is linked to the stress tensor with the transport equation (see details in Appendix C), $$\begin{array}{ll}\dot{M}(r)\hfill & =-2\pi r\Sigma {\langle u_R\rangle }_r\hfill \\ \hfill & =\frac{2\pi }{{(r^2\Omega )}^{\prime }}\left[\frac{\text{d}}{\text{d}r}\left(r^2{\displaystyle {\int }_{-H}^H{T}_{r\varphi }}\text{d}z\right)+r^2{T}_{z\varphi }{|}_{-H}^H\right]\hfill \\ \hfill & =\frac{2\pi }{{(r^2\Omega )}^{\prime }}\left[{(r^2\Sigma c_{\text{s}}^2{\alpha }_r)}^{\prime }+\frac{r^2}{2H} \Sigma c_{\text{s}}^2{\alpha }_z\right]\hfill \\ \hfill & \approx \sqrt{\frac{2\pi \Sigma r}{G}}{c}_{\text{s}}^2\left({\alpha }_r+\frac{r}{2H}{\alpha }_z\right),\hfill \end{array}$$(14)

where ′ signifies the derivative with respect to r. The approximation is obtained by making several assumptions such that d ∕d r ≈ 1∕r and Ω² ≈ G∫ 2πrΣdr ≈ 2πGΣr², which are not too far from reality (see, e.g., Fromang et al. 2013), with the purpose of allowing an appreciation of how the two terms in parentheses compare to one another. More precisely, α_r is measured in the bulk of the disk, while α_z is measuredon the surfaces, and their effects on the mass accretion rate differ by a factor r∕(2H). These two values should therefore not be compared directly.

While the stress tensor is defined along a thin circular line, the average is in practice taken within a finite volume. The α values are thus calculated as $$\begin{array}{llll}{\alpha }_{\text{Rey},r}\hfill & ={\widehat{T}}_{\text{Rey},r\varphi }/\widehat{P},\hfill & {\alpha }_{\text{Rey},z}\hfill & =\left({\widehat{T}}_{\text{Rey},z\varphi }^H-{\widehat{T}}_{\text{Rey},z\varphi }^{-H}\right)/\widehat{P},\hfill \\ {\alpha }_{\text{Max},r}\hfill & ={\widehat{T}}_{\text{Max},r\varphi }/\widehat{P},\hfill & {\alpha }_{\text{Max},z}\hfill & =\left({\widehat{T}}_{\text{Max},z\varphi }^H-{\widehat{T}}_{\text{Max},z\varphi }^{-H}\right)/\widehat{P},\hfill \\ {\alpha }_{\text{Grav},r}\hfill & ={\widehat{T}}_{\text{Grav},r\varphi }/\widehat{P},\hfill & {\alpha }_{\text{Grav},z}\hfill & =\left({\widehat{T}}_{\text{Grav},z\varphi }^H-{\widehat{T}}_{\text{Grav},z\varphi }^{-H}\right)/\widehat{P},\hfill \end{array}$$

where the hat symbol signifies the volume average. The region of the summations is an annulus of radial thickness dx, while its vertical thickness is to be specified.

The volume-averaged stress tensors and pressure are defined as $$\begin{array}{ll}{\widehat{T}}_{p\varphi ,\text{Rey}}{\displaystyle \sum \text{d}}V\hfill & ={\displaystyle \sum \delta }{u}_p\delta u_{\varphi }\rho \text{d}V={\displaystyle \sum \delta }{u}_p\delta u_{\varphi }\text{d}m\hfill \\ {\widehat{T}}_{p\varphi ,\text{Max}}{\displaystyle \sum \text{d}}V\hfill & =\frac{1}{4\pi }{\displaystyle \sum B_p}{B}_{\varphi }\text{d}V={\displaystyle \sum u_{A,p}}{u}_{A,\varphi }\text{d}m,\hfill \\ {\widehat{T}}_{p\varphi ,\text{Grav}}{\displaystyle \sum \text{d}}V\hfill & =\frac{1}{4\pi G}{\displaystyle \sum g_p}{g}_{\varphi }\text{d}V={\displaystyle \sum u_{G,p}}{u}_{G,\varphi }\text{d}m,\hfill \\ \widehat{P}{\displaystyle \sum \text{d}}V\hfill & ={\displaystyle \sum P_{\text{therm}}}\text{d}V={\displaystyle \sum c_{\text{s}}^2}\text{d}m,\hfill \end{array}$$

where $$\begin{array}{ll}{u}_A\hfill & =\frac{B}{\sqrt{4\pi \rho }},\hfill \\ u_G\hfill & =\frac{\nabla \Phi }{\sqrt{4\pi G\rho }}=\frac{-g}{\sqrt{4\pi G\rho }}.\hfill \end{array}$$

dV and dm represent the volume and mass of each cell.

The three radial terms α_Rey,r, α_Max,r, and α_Grav,r are related to the radialtransport within the disk. Their values are therefore calculated as the mean across the disk vertical extent. The summations were made within H or 3H. The total α_r is the sum of all the three components.On the other hand, α_Rey,z, α_Max,z, and α_Grav,z are the vertical counterparts that are linked to the accretion or angular momentum transport across the disk surfaces. The sum of the three terms, α_z, also leads to radial mass accretion. The summation was performed within a vertical thickness of dx around H and 3H in this case for the stress tensors, while the pressure was still averaged within the disk vertical extent. With the definition in Eq. (14), it naturally justifies the choice of normalization by the bulk pressure, whether the stress tensor is averaged in the bulk (T_rϕ) or only on the surface (T_zϕ).

Figure 14 shows the r-(top panel) and z-(lower panel) components of α. The gravitational term is mostly noise, and we do not show it for conciseness (see Appendix D for more details). The solid lines are the values measured at H and the dashed lines are those measured at 3H. For the α_r component at H, the Maxwellterm in the inner part of the disk is almost zero because the disk is strongly demagnetized as a consequence of the ambipolar diffusion (see also Sect. 5.4.2). It grows steadily when the magnetized outer part of the disk is approached and reaches ~ 10⁻³ at R_kep. The Reynolds term fluctuates around zero with a maximum absolute value lower than 10⁻³ throughout the disk. Both components grow significantly in the outer zone of the disk between R_kep and R_mag. The Reynolds component flips between positive and negative values, showing a complex transport pattern in the inner disk, which is also seen in the radial velocity (Sect. 5.4.1). The values at H and 3H are generally comparable, implying that the radial stress tensors do not cause significant stratified behavior within the vertical extent of the disk. On the other hand, the values of α_z show more significant variation across different altitudes, and this is the dominant component that leads to radial accretion within the disk. The Maxwell term at 3H is higher that at H above r ~ 5 AU, implying that most of the disk radial accretion occurs in the upper layers of the disk due to magnetic angular momentum transport. At H, the Reynolds component is significantly negative, which might imply viscous outward transport near the disk midplane.

Equation (14) implies that α_z has a stronger effect than α_r by roughly a factor 5 in the case of aspect ratio 0.1 if α_z ≃ α_r. With the slightly higher value of α_z compared toα_r, we can inferthat most of the radial mass accretion is due to the angular momentum transport along the field lines on the surface of the disk, and this accretion occurs mostly in the upper layers (between H and 3H, or even above) of the disk.

The mass accretion rate derived from α is shown in Fig. 25 with dotted curves and is compared to the directly measured values later. Strong fluctuations are present, while the overall magnitude is comparable with the directly measured values of radial accretion (see Sect. 5.5.4).

The above measurements of α along given surfaces reveal important information: The transport onto and within the disk is highly structured. To give a better idea of the transport pattern and the behavior of different stress tensors, we show α values on the r − z plane for several snapshots of the canonical run in Fig. 15. We also examine the effect of high resolution restarts in Figs. 16 and 17. Positive values are shown in red, corresponding to radial accretion (or angular momentum removal, see Eq. (14)), and negative values in blue signify radial expansion. The total effect receives contributions from all the four terms shown in the figure.

Figure 15 shows that at early times α_Rey usually dominates within the disk but does not show a regular pattern (first two snapshots). This is probably linked to the strong accretion that constantly perturbs the disk, and we recall that the disk size fluctuates in time under the effect of irregular infall of the envelope. Moreover, the disk size is very small at the beginning (10 AU for 1 AU resolution), and a good determination of the interior structure of the disk is in principle not possible. At later times, the disk grows to a few dozen AU in size, which results in a more resolved internal structure.

The middle layers of the disk do not have very well behaved values of α and some parts of the disk midplane can even expand. Most of the accretion is channeled through a high-altitude layer of the disk and moves radially inward before it reaches the surface of the disk. This is most clearly shown in the figures of α_Rey,z and α_Max,z as a dark red region slightly above the line of 3H.

Figures 16 and 17 show the comparisons with the canonical run of R_40ky_ℓ18 and R_80ky_ℓ18. At 40 kyr, the restart with increased resolution allows us to study the internal structure of the disk without having to evolve the whole simulation at high resolution since the beginning. The difference in α_Rey,r is more remarkable but is not the dominant term. Accretion occurs mostly at high altitude (see α_Max,z), and the middle layers of the disk are highly structured. At 80 kyr, the disk grows much larger, and thus even the canonical run allows satisfying evaluation of α. The similarity of the α patterns in the two runs with a difference of a factor 16 in resolution is a validation of numerical convergence. Moreover, the high-resolution restart also allows us to probe much smaller radii and shows that the similar disk and accretion structure continue at smaller radii.

Fig. 13 Plasma beta, the ratio between thermal and magnetic pressures, β = Ptherm∕Pmag, evaluated within one (orange) and three (green) times the disk scale height. The disk is highly demagnetized because ambipolar diffusion leaves the field lines outside during the collapse. The β value passes through unity around the junction between the disk and the envelope. The difference between the two curves comes from the lower magnetization level close to the disk midplane.

Fig. 14 Disk α for R_40ky_ℓ18 at 103 kyr. Top panel: Reynolds and Maxwell stresses in r direction. Bottom panel: two components in z direction. Both positive and negative values with |α| > 10−4 are plotted in logarithmic scale, and this is connected by a linear range in between.

5.4.4 Toomre Q parameter: stability and fragmentation

The Toomreparameter, Q, indicates whether the disk is subject to gravitational fragmentation. It is defined as $$Q(r)=\frac{\kappa (r)c_{\text{s}}(r)}{\pi G\Sigma (r)}\approx \frac{\Omega (r)c_{\text{s}}(r)}{\pi G\Sigma (r)},$$(24)

where $${\kappa }^2=\frac{1}{r^3}\frac{\text{d}(r^4{\Omega }^2)}{\text{d}r}=\frac{2\Omega }{r}\frac{\text{d}(r^2\Omega )}{\text{d}r}$$(25)

is the epicyclic frequency. The approximation, κ ≈ Ω, is valid for a disk in Keplerian rotation, therefore valid only within ~ R_mag. The Toomre parameter compares the stabilizing radial shear from differential rotation to its self-gravity that induces local collapse. For a value of Q > 1, the disk is stable against self-gravitating fragmentation.

Figure 18 shows the Toomre Q profile for some snapshots. The value is very high near the center because of the high temperature, and it decreases to reach the lowest point near R_kep. The value of Q increases again in the exterior magnetized part of the disk, while the approximation in Eq. (24) starts to be invalid. With the high value of Q, the whole disk is very stable against self-gravitating fragmentation. All values of c_s, Ω, and Σ are evaluated locally within H or 3H. The difference in the vertical direction mostly arises because not the entire mass is taken into account when only the region within H is considered.

The Toomre parameter Q ~ 10 near the edge of the Keplerian disk, and it reaches very high values at smaller radii (10² −10³). This extreme stability, from high temperature and/or low surface density, should be regarded with caution, and we point out some possible origins. First, the disk might be overheated with our prescription of protostellar irradiation (see Sects. 5.3 and 6). Second, the disk global profile is likely to be governed by the star-disk boundary condition, or the sink accretion algorithm from a numerical aspect. Our simulation has a percent-level disk-star mass ratio, which is somewhat lower than in some other studies (e.g., 30–40% was found by Tomida et al. 2017). This result is directly linked to the density threshold for the sink accretion (see Sect. 3.3). This has a direct consequence on the surface density of the disk, and therefore high Q values could occur if Σ is too low.

Fig. 15 α evolution of R_ℓ14 (increasing time from left to right). The four components αRey,r, αMax,r, αRey,z, αMax,z (from top to bottom). The r components are vertically averaged values between zero and a given altitude, and the z components are local values. The white curves trace the thermal scale height H and 3H. Positive values (red) can be roughly interpreted as a radial accretion, and negative values (blue) as radial expansion. The resolution is 1 AU. As the disk grows in size, the interior of the disk is described with more cells and there is a clearer pattern of α. The early disk has relatively strong αRey compared to its magnetic counter part. It is probably under constant perturbation by the strong infall from the envelope. The first two chosen snapshots are not necessarily representative of this epoch because the behavior of αRey is very irregular in time, leading to a disk that fluctuates in size around 10 AU. The infall from the envelope weakens later, and the disk becomes more regular, with its inner parts readily accreting. This accretion is mostly caused by vertical magnetic transport of angular momentum within the envelope.

5.5 Properties of the flow: infall from the envelope and accretion within the disk

The primary goal of this study is to follow the building phase of the protoplanetary disk: how mass arrives onto the disk and how mass is transported within the disk. We therefore start by considering some general properties of the global flow, and then continue with a more detailed quantification of the mass flow properties.

5.5.1 Global properties: Mass transport

Figure 19 shows R_ℓ14 (top) and R_40ky_ℓ18 (bottom) at103 kyr. Colored streamlines present the azimuthally averaged poloidal mass flux, ${\langle \rho u_p\rangle }_{\varphi }$ in addition to the specific angular momentum, ${\langle \rho r{u}_{\varphi }\rangle }_{\varphi }/{\langle \rho \rangle }_{\varphi }$, shown in gray scale. White curves mark the disk scale height H and 3H up to R_kep. The general behavior is not strongly affected by the resolution. The flow pattern varies slightly at different instants. Nonetheless, the patterns all show some similar characteristics that we present below: infall, outflow, shock, and accretion (or sometimes expansion or decretion). Figure 20 presents a schematic view of the various zones.

• Infall: The envelope is globally infalling. The infall is mostly radial near the equatorial plane. Slightly above the disk surface, the vertical mass flux increases. Most importantly, the flux intensity varies (change of color) with altitude because of the small radial component that brings the mass flux to smaller radii as it approaches the disk surface from the top. Therefore most of the infalling mass reaches small radius of within a few AU before penetrating the disk surface. We recall the source function, which is the mass flux crossing the disk surface, that we aimed to measure. This type of infalling pattern gives a surface flux that is centrally concentrated and strongly depends on the reference surface (altitude).

• Outflow: In the axial direction, a low-density outflow cavity is launched at 10–20 AU above the disk midplane, and the radial extent increases with the vertical distance from the disk. The transition between infall and outflow occurs at high altitudes above the disk, roughly near the cone z ≈ r. Immediately below this is a zone with low angular momentum (darker in the gray image), through which significant infall is channeled. There might be some circulation and recycling of the wind material across the cavity wall.

• Shock:We recall from Fig. 12 that the radial velocity in the equatorial plane reaches a fraction of the free-fall velocity and quickly drops at ~R_mag. In the post-shock region, the rotation replaces the infall as the dominant kinetic energy. A slowly expanding zone appears very often in the magnetized outer region of the disk and meets the radial infall at the shock front. This shock is located at r ~ 15 AU in both panels of Fig. 19, where the radial velocity apparently decreases. In R_ℓ14, it is more evident that the gas is pushed upward due to the increased pressure and falls back to the disk surface at smaller radius, while the streamlines in R_40ky_ℓ18 are more regular and there is only very mild upward motion when the gas reaches the shock front. The radial flow inside the disk is sometimes outward, causing circulation at the edge of the disk: matter moves upward at the shock and falls back from the top. This behavior has been suggested in observations (Sakai et al. 2017), and the shock might be traced with chemical molecules. However, in our case, the disk is probably too thin and the shock-generated heat is quickly radiated away.

• Accretion: Within the disk (z ≲ 3H), the behavior of the flow is less clear. R_ℓ14 shows an almost vertical flux downward, while in R_40ky_ℓ18 the flux is upward. The disk internal dynamics does not converge numerically, and conclusions cannot be drawn. Globally speaking, the innermost parts of the disk (a few AU) are always in radial accretion, which leads to the mass growth of the central star. However, the radial velocity within the rest of the disk is not strictly negative all the time. Outward expansion sometimes occurs in the middle layers of the disk. As mentioned earlier, this occurs very often in the magnetized zone of the disk (between R_kep and R_mag). However, the disk is confined by the magnetized zone and the infalling envelope and thus cannot expand freely.

At 103 kyr, the disk radius is roughly 15 AU, while it grows to 50 at 138 kyr (Fig. 21). The general behavior is the same as described above, and the region with low angular momentum becomes a clear channel of mass infall. The expanding magnetized zone is very large and almost comparable to the disk radius, R_kep. Here again, the disk internal dynamics is not identical in the two runs with different resolutions. Nonetheless, the two main remarks remain valid: matter might circulate near the wind-infall wall and at the disk edge.

Fig. 16 Measured α values of R_ℓ14 (left) and R_40ky_ℓ18 (right) at 103 kyr. The disk radius is about 15 AU. The bins for averaging are chosen to be larger than dx in the high-resolution run in order to remove numerical noise. These figures clearly show that the technique of restart using higher resolutions is needed to study the interior properties of the disk. The z components are more important because they dominate in governing the mass accretion inside the disk. The middle layers of the disk are more irregular with some partially expanding parts, while most of the accretion occurs through the upper layer slightly above 3H and reaches a small radius before penetrating the surface of the disk.

Fig. 17 Measured α values of R_ℓ14 (left) and R_80ky_ℓ18 (right) at 138 kyr. Disk radius (Rkep) grows to 50 AU at this time. The values are evaluated on a grid of 1 AU resolution for both runs. The global patterns are compatible for the two runs with 16 times difference in resolution. The middle layers of the disk are radially structured and partially expending at times, while the accretion occurs through a channel slightly above 3H and brings a significant amount of mass from the infalling envelope to a small radius without passing through the disk.

5.5.2 Global properties: angular momentum transport

The commonly accepted scenario describes the loss of angular momentum due to various transport mechanisms in the envelope, which allows the matter to be accreted within the disk. The equation of angular momentum transport is obtained by multiplying the ϕ component of the momentum conservationequation by r and averaging azimuthally, which gives $$\frac{\partial }{\partial t}{\langle \rho r{u}_{\varphi }\rangle }_{\varphi }+\nabla \cdot \left[{\langle r\rho u_p{u}_{\varphi }-r\frac{B_p{B}_{\varphi }}{4\pi }\rangle }_{\varphi }\right]=0.$$(26)

The second term is the divergence of the poloidal flux of angular momentum. By separating the toroidal velocity into the circular and fluctuating components u_ϕ = rΩ + δu_ϕ, where $r\Omega ={\langle \rho u_{\varphi }\rangle }_{\varphi }/{\langle \rho \rangle }_{\varphi }$, we can rewrite the angular momentum flux as $$\begin{array}{ll}{\mathcal{F}}_j\hfill & ={\langle r\rho u_p{u}_{\varphi }-r\frac{B_p{B}_{\varphi }}{4\pi }\rangle }_{\varphi }\hfill \\ \hfill & =r^2\Omega {\langle \rho u_p\rangle }_{\varphi }+r{\langle \rho \delta u_p\delta u_{\varphi }\rangle }_{\varphi }-r\frac{{\langle B_p{B}_{\varphi }\rangle }_{\varphi }}{4\pi }.\hfill \end{array}$$(27)

The first term is the laminar transport, which describes the angular momentum flux associated with the mass flux if the momentum is to be frozen within the mass. The second term is the turbulent transport term that originated from the cross correlation of the fluctuations. This is the main transport mechanism in the α-disk model that allows the disk to accrete. The last term is the magnetic transports associated with the Maxwell stress tensor, which tends to deviate the movement along the field line direction, thus accelerating or braking the rotation.

In the case of an isolated disk at later phases, the focus lies more on the angular momentum within the disk, and there is no important net mass flux. The tensor terms are therefore of main interest. In this study, however, the disk actively receives mass from the infalling envelope and the laminar transport must not be neglected. Figure 22 shows the laminar, turbulent, magnetic, and total angular momentum fluxes (from top to bottom) at 103 kyr at two resolutions (left and right). From the pattern of angular momentum transport, the disk-envelope system can be divided into four zones (see Fig. 20 for an illustration).

• The braked infalling layer: this region is dominated by the vertical upward transport of angular momentum by the magnetic field while the velocity is downward. This means that during the infall, the gas is braked by the magnetic field that is not aligned with the flow. This is the dark zone in Fig. 19 where the specific angular momentum is lowest and the mass is strongly infalling.

• The freeinfalling layer: the region contained within a few scale heights above the disk surface is dominated by the free infall of mass. The turbulent and magnetic stress tensors transport the angular momentum upward, while this is dominated by the laminar flow that directly advects the angular momentum associated with the gas. In this region, the magnetic field is mostly lost through ambipolar diffusion. This zone is clearly seen by comparing Figs. 19 and 22, where the mass flux and the total angular momentum flux are almost in the same direction.

• The wind cavity: this low-density region is dominated by the vertical upward transport of angular momentum by the magnetic field. Matter moves outward along field lines, and its angular momentum increases.

• The disk interior: the disk interior (below a few H) is dominated by the laminar transport of the angular momentum. This is linked to the global behavior of the disk and has to be investigated with great care. While radially outward magnetic transport and vertical Reynolds transport exist as well, these two contributions are weak. We refrain from overinterpreting the results because the low- and high-resolution runs are not consistent at this scale.

Throughout thebraked layer and the wind cavity, the magnetic angular momentum flux is almost constant. This means that the field lines channel out the angular momentum in an almost stationary way and eventually launch it away through the wind.

Figure 23 shows the same plots at 138 kyr, where the disk grows to a larger size. The generalbehavior is the same, while we have better resolved disk internal dynamics, especially in R_80ky_ℓ18. The disk midplane at this moment expands. It is clear that the angular momentum transport is largely dominated by the laminar transport that simply advects the angular momentum along with the mass. These results show that the α-disk model is probably not suitable for an embedded disk because the disk is not in a stationary state and is very sensitive to irregularities in the infalling envelope.

Fig. 18 Toomre Q ≈ Ωcs∕(πGΣ) parameter evaluated within one (orange) and three (green) times the disk scale heights. The same snapshots are taken as in Fig. 4. The disk is very stable against self-gravitating fragmentation. The difference between the two curves mainly arises because the total mass is not taken into account for the value within H and 3H.

Fig. 19 Poloidal mass flux (${\langle \rho u_p\rangle }_{\varphi }$) and specific angular momentum ($r{\langle u_{\varphi }\rangle }_{\varphi }$) around the disk region (color streamlines in the gray image, in SI units), for R_ℓ14 (top) and R_40ky_ℓ18 (bottom)at simulation time 103 kyr. All quantities are azimuthally averaged and averaged between the upper and lower planes. White curves outline the disk scale height H and 3H at r < Rkep. An outflow is launched at high altitudes at z ≈ 20 AU in R_ℓ14. In the high-resolution run, the flow varies slightly more, but it is qualitatively similar. The outflow is launched closer to the disk (z ≈ 10 AU) and the disk size is slightly smaller. The outflow has relatively low density and manifestsitself as a cavity. Significant infall from the envelope is channeled through a zoneof low angular momentum at several scale heights (darker gray region near z ≈ r) of the disk. The horizontal infall from the envelope reaches the outer edge of the disk (shown in Fig. 2 at r ≈ 15 AU), and thedensity contrast generates a shock that pushes the gas to higher altitude. This creates a counterclockwise poloidal circulation near the disk edge between 10 and 15 AU in R_40ky_ℓ18 because the outer disk decretes at this instant, while it is not seen in the upper panel.

Fig. 20 Zones showing the main properties of the disk-envelope system in Sect. 5.5.1 and 5.5.2.

Fig. 21 Poloidal mass flux (${\langle \rho u_p\rangle }_{\varphi }$) and specific angular momentum ($r{\langle u_{\varphi }\rangle }_{\varphi }$) around the disk region (color streamlines in the gray image, in SI units) for R_ℓ14 (top) and R_80ky_ℓ18 (bottom) at simulation time 138 kyr. All quantities are azimuthally averaged and averaged between the upper and lower planes. The white curves outline the disk scale height H and 3H at r < Rkep. An outflow is launched at z ≈ 20 AU. The disk has grown in size, while the flow pattern stays very similar to that in Fig. 19. At this moment, the magnetized zone is almost twice the size of the Keplerian disk (Rmag ≈ 2Rkep), and the counterclockwise poloidal circulation around Rkep is evident inboth panels and occupies almost the entire magnetized zone and the outer Keplerian disk.

Fig. 22 Angular momentum flux and angular momentum ($r{\langle \rho u_{\varphi }\rangle }_{\varphi }$) around the disk region (color streamlines in the gray image, in SI units), for R_ℓ14 (left) and R_40ky_ℓ18 (right) at 103 kyr. All quantities are azimuthally averaged and averaged between the upper and lower planes. From top to bottom: Laminar transport ($r{\langle u_{\varphi }\rangle }_{\varphi }{\langle \rho u_p\rangle }_{\varphi }$), turbulent transport ($r{\langle \rho \delta u_{\varphi }\delta u_p\rangle }_{\varphi }$), magnetic transport ($-r{\langle B_{\varphi }{B}_p\rangle }_{\varphi }/(4\pi )$), and their sum. White curves outline the disk scale height H. The laminar flow convects mass with positive angular momentum, while the tensor terms (Reynolds and Maxwell) transport the angular momentum outward. In particular, the magnetic tensor dominates in the envelope, and this allows the mass to fall onto the disk through the loss of angular momentum. The general pattern persists with improved resolution. White patches in the gray image are due to locally negative angular momentum. This is an artifact of the averaging because the selected axis is not exactly aligned with the asymmetric outflow.

Fig. 23 Angular momentum flux and angular momentum ($r{\langle \rho u_{\varphi }\rangle }_{\varphi }$) in the disk surrounding (color streamlines in the gray image, in SI units) for R_ℓ14 (left) and R_80ky_ℓ18 (right) at 138 kyr. The disk scale height H is shown in white. All quantities are azimuthally averaged and averaged between the upper and lower planes. From top to bottom: Laminar transport ($r{\langle u_{\varphi }\rangle }_{\varphi }{\langle \rho u_p\rangle }_{\varphi }$), turbulent transport ($r{\langle \rho \delta u_{\varphi }\delta u_p\rangle }_{\varphi }$), magnetic transport ($-r{\langle B_{\varphi }{B}_p\rangle }_{\varphi }/(4\pi )$), and their sum.

5.5.3 Infall onto the disk: source function

The source function, S(r), in units of mass per unit time per unit surface, describes the rate and the spatial distribution of mass arrival from the envelope onto the disk surface. The surface we refer to when we discuss the mass flux is always the projection of the disk surface onto the equatorial plane, which is equal to the actual (inclined) disk surface multiplied by the cos (θ) of the local inclination angle. This choice is made such that the description is consistent with the classical model that regards the disk as infinitely thin (e.g., Hueso & Guillot 2005).

The disk scale height h_p (see Sect. 5.1.1) was used to define the surface at which the source function was evaluated. To avoid numerical errors due to the small-scale fluctuations of h_p, we first fit the scale height to a power-law expression, $$H=\eta {\left(\frac{r}{1 \text{AU}}\right)}^p\text{}\text{},$$(28)

such that it became a monotonically increasing smooth function. The fitted parameters are shown in Fig. 4. The disk aspect ratio is almost constant at 10%. We chose to measure the source function at h_s = H and h_s = 3H. The source function is therefore measured as $$\begin{array}{ll}S(r)\hfill & =-\rho (r,z)\left[u_z(r,z)-u_r(r,z)\tan (\theta (r))\right]\hfill \\ \hfill & =-\rho (r,h_{\text{s}}(r))\left[u_z(r,h_{\text{s}}(r))-u_r(r,h_{\text{s}}(r))\frac{\text{d}{h}_{\text{s}}}{\text{d}r}(r)\right].\hfill \end{array}$$(29)

Figure 24 shows snapshots of the measured S(r) at H (orange) and 3H (green) at two resolutions. Closer to the equatorial plane, the source function is more centrally concentrated, as we discussed earlier; this is due to the converging flow pattern (see Sects. 5.4.3 and 5.5.2). When we compare this with the source functions from Ulrich (1976) and Hueso & Guillot (2005), the measured S(r) is much steeper, meaning that the infall onto the disk is more centrally concentrated for our simulated disk. Particularly, this means that most of the disk at r ≳ 3−4 AU does not receive much material and a significant amount of mass falls straight onto the central star. We note that the choice of reference plane for the mass flux measurement is not a simple matter. The lower panel of Fig. 24 shows that the mass flow crosses the disk surface and moves upward, giving negative flux. A reference surface therefore has to be sufficiently far away from the midplane to avoid being affected by the disk internal motions.

Moreover, because the size of the sink particle accretion radius is finite, the flux within the disk scale height at small radius is not physical. We therefore measured the flux at the horizontal surface z = 4dx (shown in gray). It is evident that the extreme high flux values at small radii are not physical. Very close to the disk center, the measurement of the vertical and horizontal flux becomes sensitive to the exact selection of the axis, and there might be some confusion between the two. When we calculated the total mass flux, we therefore used min (4dx, H) for the disk surface. This introduces some uncertainties in the measurement, but it is the best we can do at this given resolution with the constraints from the algorithm.

5.5.4 Accretion inside the disk

The straightforward measure of the mass accretion rate, Ṁ_r, is the mass that crosses radius r inward per unit time. This is calculated as $${\dot{M}}_r(r)=-2\pi r{\sum }_{-h_{\text{s}}}^{h_{\text{s}}}{\langle \rho u_r(r,z)\rangle }_{\varphi }\text{d}z,$$(30)

where ρu_r is the local radial mass flux. This radial accretion rate is shown in Fig. 25 with dashed curves for the integrated values within H (orange) and 3H (green). At small radii, the radial flux within ±4dx, which is larger than the scale height, is shown, and therefore the green curves are superimposed on the orange curves. As discussed earlier, the disk can sometimes be in partial expansion (negative Ṁ_r (r)), while in general, the inner part of the disk (below ~5 AU) is normally accreting, with small local fluctuations. The discrepancy between the H and 3H curves shows that the radial accretion within the disk mostly occurs in the upper layers. This is shown in the two lower panels of Fig. 25. Beyond ~ 10 AU, the green (3H) curve is often above the orange (H) curve, meaning more accretion (positive) or less decretion (negative). In the region slightly beyond the sink accretion radius (second vertical gray line), the value of Ṁ_r(r) increases toward the center of the disk in all four panels, meaning that either the disk is being emptied or that there is a surface source.

To examine this, we also calculated the surface mass accretion rate, that is, the source function integrated over the disk surface. We chose to perform this integration between 0 and r, giving $${\dot{M}}_s(r)={\sum }_0^{r}S(r)2\pi r\text{d}r.$$(31)

This is shown in the same figure with solid curves for source functions measured at H (orange) and 3H (green). Because of the numerical limits described above, we used the flux at the surface z = 4dx when H was smaller than this value. For R_ℓ14 at 103 kyr, the disk is so small that H < 4dx for almost the entire disk.

Most of the mass penetrates the disk surface where the slope of Ṁ_s (r) is the steepest. The total mass that falls onto the disk, with significant contribution from this zone at a few AU, is broadly comparable with the mass that falls radially inward onto the star (Ṁ_r(r) at small radii).As a general remark, the infalling material from the envelope can reach quite a small radius of a few AU before it penetrates the disk surface. Inside the disk, the inward accretion is well established only for the innermost of the disk, while some of the mass that falls from the surface onto the disk might move outward. This expansion reaches a shock generated by the radial infall form the envelope near the equatorial plane and brings the material to high altitudes, possibly generating a circulation in the outer part of the disk.

Radial accretion behavior is present in all panels in the marginally resolved region, that is, a few times the sink accretion radius. Nonetheless, we would like to caution that this zone is primarily dominated by the numerical viscosity, and this is probably why the accretion is readily established. This should not have strong consequences on the long-term stellar mass growth because the infall from large scales is the dominant factor. However, this numerical effect does determine whether the evolution of the star-disk system can be described on short timescales, which remains a computational challenge (e.g., Kuffmeier et al. 2018). In the outer part of the high-resolution runs, both restarts show oscillating radial motions. This again supports our proposition that the disk is highly sensitive to the environment, while the role of the disk as a mass buffer between the envelope and the star is less evident.

Figure 25 also shows Ṁ_α, the mass accretion rate estimated from the α parameter (see Sect. 5.4.3, Eq. (14)). The values of Ṁ_r and Ṁ_α do not match perfectly. However, they are similar in magnitude.

Fig. 24 Source function S(r) measured at H (orange) and3H (green) in R_ℓ14 and R_40ky_ℓ18 at 103 ky, R_ℓ14 and R_80ky_ℓ18 at 138 ky (from top to bottom). The flux measured at the horizontal plane z = 4d x is shown in gray. The curves in lighter colors are the analytical functions from pure hydrodynamic predictions that have the same amount of total flux (see Fig. 3). The source function measured in the simulations is much steeper than that of the classical model, and strong infall occurs in the central part of the disk. At a given radius, the flux is generally greater at lower altitudes, indicating that there is a strong radial flux in the interior of the disk. At small radii (below the sink accretion radius), the measured flux at H∕3H is no longer realistic, and we used the measurement at z = 4dx to calculate the total flux.

Fig. 25 Massaccretion rate measured across the disk surface (Ṁs(< r)), measured radially inside the disk (Ṁr(r)), and inferred from the measured α (Ṁα (r)) for R_ℓ14 and R_40ky_ℓ18 at 103 ky, R_ℓ14 and R_80ky_ℓ18 at 138 ky (from top to bottom). Both positive and negative values with |Ṁ | > 10−7 are plotted in logarithmic scale, and this is connected with a linear range in between. Most of the surface flux arrives where the slope of Ṁs(< r) is the steepest. Below this radius, the disk readily accretes in the radial direction only at a very small radius (<1 AU), and the outer part can have complex motions. The radial mass flux derived from α is not exactly the same as Ṁr(r), although the amplitudes are comparable.

6 Comparison to standard alpha-disk models

The present paper shows that the disk formed during the infall of the dense core has a number of similarities to but also differences with classical isolated disks assumed in standard alpha-disk models. We report the most notable similarities below.

• 1.

  The Keplerian disk structure is established rapidly in less than 61 kyr, a timescale shorter than the envelope infall timescale.

• 2.

  The disk mass is about 1% of the stellar mass after about 80 kyr (Fig. 7). This is a typical disk-to-star mass ratio (like for the MMSN, see, e.g., Hayashi 1981). However, we note that the stellar mass is only 0.2 M_⊙ at the end of the simulation so that we cannot confirm at the moment that this results is valid up to one solar mass. It is therefore premature to state that there is enough mass to build the S.S. Moreover, we caution again that the disk mass is (unfortunately) sensitive to the disk inner boundary condition (n_acc).

• 3.

  The midplane of the disk is almost devoid of MHD turbulence with α_r < 10⁻⁴. This is in agreement with previous studies arguing that the midplane of a protoplanetary disk should be dead because of lack of ionization (see Turner et al. 2014, for a review). Several authors showed that low values of α in the disk midplane can lead to efficient planetesimal formation (Drążkowska & Dullemond 2014, 2018; Charnoz et al. 2019). The resulting disk might be the birthplace of a massive planetesimal population. However, we caution that the Toomre Q parameter is high in our simulation. Planetesimal formation can be allowed if Σ is higher in reality (see Sect. 5.4.4), or if the metallicity is enhanced toward the disk midplane with dust enrichment as high as a factor ≃2, as indeed suggested recently by Lebreuilly et al. (2019, 2020).

We report the notable differences with classical alpha-disk models below.

• 1.

  The disk in our simulation is truncated at about 20 AU through magnetic braking, consistent with Hennebelle et al. (2016). This may imply that the classical picture of a viscously spreading disk up to several dozen or 100 AU is maybe not the general rule, whereas a disk like this does exist (see, e.g., Isella et al. 2010). However, recent ALMA observations showed that the majority of class 0 disk are likely small (Andrews et al. 2018; Maury et al. 2019). The absence of viscous spreading may preclude the outward transport of CAIs because of outward viscuous relaxation of the disk (Yang & Ciesla 2012; Pignatale et al. 2018). However, alternative solutions could exist because transient outward-directed flows are detected in the disk (see Sect. 5.5) and might drive CAIs outward from the protostar. Alternatively, the disk might expand in a later phase when the envelope is almost exhausted.

  Interestingly, the truncation of the disk at 20 AU agrees with some versions of the Nice model, as well as the more generalized giant planet instability model that appeared later. In the recent view of S.S. formation, the giant planets formed in a much more compact configuration than today, within about 15–20 AU of the Sun (see, e.g., Gomes et al. 2005; Tsiganis et al. 2005; Morbidelli & Nesvorný 2020). The main argument for this is the necessity of forming a low-mass Kuiper belt with a correct orbital architecture in which an accumulation of bodies is observed in mean-motion resonance with Neptune. To reach the observed state, it is necessary to invoke an outward migration of Uranus and Neptune starting from 15–20 AU. When the gas disperses, the four giant planets are still in a compact configuration, with Jupiter and Saturn close to their 3:2 mean motion resonance. This scenario agrees with disk and planet evolutionary models (Bitsch et al. 2015), at least for Jupiter and Saturn.

  The Nice model originally aimed to explain the lunar bombardment (Gomes et al. 2005). The key idea is that the orbits of giant planets were unstable during the S.S. evolution. If the giant planets were formed in a more compact configuration than today (as suggested in the Nice model, or as a result of evolution within the protoplanetary disk, see Kley & Crida 2008; Walsh et al. 2011), then they should at some moment have experienced a global dynamical instability, which led to a very short and violent reorganization of their orbits (Deienno et al. 2017; Clement et al. 2018). During this rearrangement, small-body belts were destabilized, resulting in massive bombardment of nearby planets. For example, the outward migration of Neptune may have pushed the original Kuiper belt outward, which might initially have been contained within about 20 AU and not at the current position of about 50 AU. It is therefore possible that the young S.S. was more compact than today, and reached its current dimensions after such an instability (Levison et al. 2008).

  In the original version of the Nice model, the giant-planets instability occurred at 800 Myr after CAIs. More recent models of giant planet instability, using the configuration of the four giant planets at the end of the solar nebula obtained from hydrodynamical simulations of planet migration, also converge to a formation of Neptune within 25 AU (Nesvorný & Morbidelli 2012). de Sousa Ribeiro et al. (2020) suggested that this event could have occurred as early as in the first million years. Forming Uranus and Neptune at 30 AU and beyond would have disrupted the Kuiper belt and in particular the population of dynamically cold objects (low orbital eccentricity and inclination). This suggests that the disk itself was initially in a compact configuration and presumably did not extend to 100 AU.

  Nonetheless, the disk in our simulation also experienced a rapid grow in size after ~ 120 kyr. This suggest that the disk size is sensitive to the properties of the infalling material, which can be highly variable in time (see also Vorobyov et al. 2015; Kuffmeier et al. 2017). Therefore more constraints are needed to reconstruct the history of the S.S.

• 2.

  The formed disk is significantly hotter than classical isolated disks. The snow line starts at about 10 AU and shifts inward down to 5 AU, where it then remains constant (until the end of the simulation) while the accretion drops to about 10⁻⁷ M_⊙ yr⁻¹. This temperature profile is notably higher than in passive radiatively heated disks (see, e.g., Chiang & Goldreich 1997), where the snow line falls around 2 AU (independently of the accretion rate, and assuming a 2.5 R_⊙ star with T = 4000 K), or isolated disk with both viscous and radiative heating: for example, Bitsch et al. (2015) reported that the snow line lay between 2 and 5 AU for an accretion rate lower than or equal to 10⁻⁷ M_⊙ yr⁻¹, consistent with the detailed model of D’Alessio et al. (1998), who reported a snow line around 2 AU for an accretion rate of ~ 10⁻⁷ M_⊙ yr⁻¹. A snow line located at 5 AU would be the ideal position to trigger fast planetesimal formation of the Jupiter core, as it has been demonstrated that streaming instability and planetesimals formation is strongly enhanced at the snow line (Drążkowska & Alibert 2017; Charnoz et al. 2019).

  Nonetheless, the temperature of the disk should be interpreted with caution because it is subject to some model uncertainties. First, the prescription by Kuiper & Yorke (2013) was used to describe the protostellar irradiation on the Hosowaka track. This heating might be excessive and overheats the surrounding of the protostar. Figure 11 shows that the snow line is shifted to 2–3 AU in the absence of this heating term. Second, the radiative transfer was treated with a flux-limited diffusion scheme, which is mostly valid for infrared radiation, while optical photons may escape more easily. In reality, the absorption of optical radiation by dust and the re-emission in infrared wavelengths is not complete and needs to be further investigated with care.

• 3.

  The disk surface density does not follow the classical −1.5 exponent of the MMSN (Hayashi 1981) but is closer to −0.5, with a sharpdrop beyond R_kep. This recalls observed protoplanetary disks, see, for instance, Isella et al. (2010), who suggested that the surface densities have an exponent >−1 and that the outer edge of the disks is between 23 and 27 AU.

7 Comparison with similar numerical works

We can draw essential conclusions from previous numerical studies and our current work. Primarily, the disk sensitively responds to its environment, that is, the collapsing prestellar core. Therefore class 0/I disks should not be studied in the same way as isolated disks. Second, whether the disk has significant effects on the evolution of the star remains debatable. Bursty mass accretion may result from fragmentation within the disk, while the numerical origin of these gravitational instabilities is not yet completely ruled out. An alternative way of regarding the disk as an accessory of the protostar and an inevitable existence due to residual angular momentum is emerging. Evolution of the star is determined mostly from the core scale and is partly mediated by the disk. This deviates from the widely accepted view that mass first falls onto the disk and then the star accretes from the disk.

Vorobyov & Basu (2010, 2015) and Vorobyov et al. (2015) employed a 2D thin-disk approximation. These models cannot describe the mass accretion that reaches the star through upper layers and everything seems to transit through the disk. In the embedded phase, the stability conditions in the disk are mostly reflected in the fluctuations (bursts) of the mass accretion rate, while the accretion rates from the envelope onto the disk and that from the disk onto the star are broadly consistent (for class 0/I, where envelope infall is significant).

Kuffmeier et al. (2018) resimulated small regions around six sink particles selected from a giant molecular cloud simulation at increased resolution of 2 AU for 200 kyr. One of these runs was further evolved for 1000 yr at 0.06 AU resolution. These parameters are similar to ours, and we share some common findings. First, the overall mass accretion history is likely determined by large-scale flow properties (see also Kuffmeier et al. 2017) instead of the viscous accretion process in the disk. Second, the sink mass accretion mediated through the disk is bursty, although it might be of numerical origin. Finally, significant mass infall arrives onto the star without transiting radially through the disk.

8 Conclusions

We have simulated the collapse of a magnetized prestellar core, considering the nonideal MHD effects of the ambipolar diffusion that removes the magnetic field from the gas. The temperature was calculated with a flux-limited diffusion scheme for the radiative transfer. We followed the evolution during the first 100 kyr after protoplanetary disk formation, and we performed high-resolution restarts in parallel to study the numerical convergence as well as the detailed interior structure of the disk. This is one of the first works to try to self-consistently form a protoplanetary disk and resolve the interior structure at the same time. Nonetheless, we recall that the numerical convergence is not strictly demonstrated in this work and challenges remain.

The collapse problem is known to be computationally challenging because the dynamical ranges are wide. For this reason, the formation of the protoplanetary disk has rarely been followed for a long duration. On the other hand, the detailed structure of the disk interior is usually studied with a prescribed disk or a shearing box that contains a small region of the disk. In this study, we have followed the evolution of the disk during a relatively long period of time at 1 AU resolution. The disk is formed self-consistently from a prestellar core that is a few thousand AU in radius. This compromise in resolution allowed us to reach more evolved stages within a reasonable computation time. In order to study the detailed interior dynamics of the disk, we selected snapshots from the canonical simulation and reran with increased resolution down to 0.06 AU.

From this work, we were able to follow the time evolution of disk global properties such as its mass and size. The disk is established roughly 60 kyr after the beginning of the simulation and grows slowly in size. The mass accretion rate is a decreasing function of time, and the disk gradually approaches an isolated state. The global appearance (density structure) of the disk is quasi-stationary and has mass 1.5% of the stellar mass, which is massive enough with respect to the star to form planets. However, the dynamics is likely to be predominantly governed by infall from the envelope, which inherits sometimes irregularities from the prestellar core. As a result, the transport of matter and angular momentum is mostly laminar and regulated by the magnetized infall. The α−disk model that describes the turbulent transport, linked to the local surface density and temperature, is probably not applicable for a disk in the embedded phase.

Because of the small size of the disk (<25 AU) and the strong MHD turbulence in the outer regions, all planets must form within ~ 30 AU (R_kep). This is qualitatively consistent with the modern view of S.S. formation, where all planets are suggested to have formed well inside 20 AU (e.g., Raymond & Izidoro 2017; Morbidelli et al. 2018) and most planetesimals may have accreted up to 30 AU (Gomes et al. 2004; Nesvorny et al. 2013; Nesvorný et al. 2018). The temperature profile evolution can also be extracted and allows determining the migration of the snow line. The radial temperature profile is consistent with a power-law exponent between −1 and −0.5. This is somewhere between a purely irradiated disk and a viscously heated disk. However, the physical origin of this profile is likely the radiative diffusion in the disk radial direction instead of the aforementioned mechanisms. The snow line starts at about 10 AU and ends at about 5 AU after about 100 kyr evolution. This is much farther away than commonly assumed for isolated and passively irradiated disks. This result should be regarded with care because it might have some numerical origins. First of all, the flux-limited diffusion scheme method treats only infrared radiation and thus could trap more energy than it should in reality. Second, the temperature is highly sensitive to the density (optical depth) of the disk, which in turn depends on the inner boundary conditions of the disk. Moreover, we recall that the stellar mass has only reached ~ 0.3 M_⊙ in this study. Initial conditions with higher mass should be explored further.

The disk has an inner part with a shallow power-law surface density profile (≳ −1 exponent) and is almost in Keplerian rotation. While the inner part is warm and almost demagnetized, an outer magnetized zone links the disk to the infalling envelope. The size of this zone varies significantly in time, from a few AU to almost comparable to the radius of the disk. The pattern of mass infall from the envelope is much more complicated than that described by the classical models that conserve angular momentum. The infall arrives mostly from a channel at altitude ≳ 3H, reaching the interior of the disk without penetrating the disk surface at large radius, which was also suggested by Kuffmeier et al. (2018). The midplane of the disk, on the other hand, can sometimes expand. The expanding material from the disk midplane sometimes exceeds the Keplerian radius, R_kep, of the disk and meets the infalling material in the magnetized zone between R_kep and R_mag. Some material is pushed up vertically and falls down again through the disk surface at a smaller radius. This might generate some circulation of material, which remains to be verified using tracer particles or passive scalar fields. In particular, the source function evaluated from our simulations is much more centrally peaked than in classical models, and most of the mass that finally rests inside the star does not transit through the whole disk.

Acknowledgements

We thank the referee, Enrique Vazquez-Semadeni, for thorough comments that have significantly improved the manuscript. This work was granted access to HPC resources of CINES under the allocation x2014047023 made by GENCI (Grand Equipement National de Calcul Intensif). Y.-N. Lee acknowledges the financial support of the UnivEarthS Labex program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02), MOST (108-2636-M-003-001, 109-2636-M-003-001), and the MOE Yushan Young Scholar program.

Appendix A Identification of the disk

Because we are interested in the interior properties of the disk as well as in its accretion onto the central star, we chose the sink particle as the center of the star-disk system. Despite some non-axisymmetric oscillation of the disk, the gravitational effects from the central star cancel out when calculating azimuthally averaged quantities, α in particular. This allows a discussion of the disk properties in coherence with the axisymmetric disk models (which is the case for almost all existing models). The disk is identified in two steps. The first step is based on pure density thresholding, and the second step is morerefined and takes the radial and vertical mass distributions of the disk into account (see Sect. 5.1).

The disk is significantly marked by a density jump, which is used as a simple first identification. The rotating dense disk is surrounded by a collapsing diffuse envelope. The procedure is as follows:

• 1.

  Calculate the inverse cumulative density function (CDF) of mass against density, that is, the amount of mass that is at a density above a certain value, $$\mathcal{M}(\rho )={\displaystyle {\int }_{\rho }^{\infty }\text{d}}m.$$(A.1)

  This results in a decreasing function of ρ that has a plateau, which corresponds to the density jump between the disk and the envelope.

• 2.

  Define the disk threshold ρ_th density as the point where the derivative of the CDF is highest, that is, where the CDF is flattest. This writes $${\rho }_{\text{th}}=\arg \max \text{d}\log \mathcal{M}(\rho )/\text{d}\log \rho .$$(A.2)

• 3.

  Select all the cells that have density ρ ≥ ρ_th.

• 4.

  Calculate the angular momentum of these cells with respect to the sink particle and define this as the disk axis.

A simple density thresholding is robust to the turbulent fluctuations for most of the time, but not always. As the high-density regions dominate the mass inside the disk, the rotation axis we defined is robust to slightly different definitions of the disk. The rotation in the initial condition is aligned to the − z axis, and the axis of the disk is very close to this direction.

Appendix B Disk vertical temperature distribution

Figure B.1 shows and example of disk vertical temperature distribution, where the temperature from each cell is plotted directly. Within a few times the scale height, the disk is basically isothermal, and decrease in temperatureis only significant far away from the disk midplane. Vertical temperature variation is more clearly seen at small radii, that is, within 7 AU, where the vertical structure of the disk is not well resolved, however. The temperature dispersion among cells at same radial distance is sometimes more prominent than the vertical variation.

Fig. B.1 Temperature plotted against altitude, z, at various distance to the central star, r, from R_ℓ14 (the disk aspect ratio is ~0.1 (Sect. 5.1.1). Every dot corresponds to one cell in the simulation, without any averaging. The discrete distribution comes from the slight difference between the grid direction and the selected disk axis.

Appendix C Stress tensors linked to the transport

Here we self-coherently define the stress tensors used in this work to avoid confusion. Following Balbus & Hawley (1998), the angular momentum conservation averaged in azimuthal direction becomes $$\frac{\partial }{\partial t}{\langle \rho r^2\Omega \rangle }_{\varphi }=\nabla \cdot r[{\langle \rho r\Omega u_p\rangle }_{\varphi }+T]=0,$$(C.1)

where T is the poloidal stress tensor, $$T={\langle \rho u_{\varphi }{u}_p-B_{\varphi }{B}_p/4\pi \rangle }_{\varphi },$$(C.2)

and W = $T/{\langle \rho \rangle }_{\varphi }$. Instead of integrating over the whole vertical domain, we integrated over the disk scale height [−H, H], and defined a mass-weighted mean $${\langle X\rangle }_{\rho }=\frac{1}{2\pi {\Sigma }_H\Delta r}{\displaystyle {\int }_{-H}^H{\displaystyle {\int }_{r-\Delta r/2}^{r+\Delta r/2}{\displaystyle {\int }_0^{2\pi }\rho }}}X\text{d}\varphi \text{d}r\text{d}z,$$(C.3)

where Δr is a radial averaging scale chosen such that the profile is not significantly affected by fluctuations. The choice of Δr = dx is reasonable for most variables, while for the stress tensors, higher values are used, and this is discussed in Appendix D. It is easily shown from this definition that the disk surface density within a the thickness ± H, Σ_H, can be found by calculating ${\langle 1\rangle }_{\rho }=1$.

Based on this averaging, the equation of mass conservation becomes $$\frac{\partial {\Sigma }_H}{\partial t}+\frac{1}{r}\frac{\partial (r{\Sigma }_H{\langle u_r\rangle }_{\rho })}{\partial r}+{\langle \rho u_z\rangle }_{\varphi }^H-{\langle \rho u_z\rangle }_{\varphi }^{-H}=0,$$(C.4)

where the superscript ±H denotes the values at ±H. The second term is the radial mass flux, and the last two terms are the flux across the disk surfaces.

Likewise, Eq. (C.1) becomes $$\begin{array}{l}\frac{\partial }{\partial t}\left({\Sigma }_H{r}^2{\langle \Omega \rangle }_{\rho }\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r^3{\Sigma }_H{\langle \Omega u_r\rangle }_{\rho }+r^2{\Sigma }_H{\langle W_{r\varphi }\rangle }_H\right)\\ +r^2\left({\Omega }^H{\langle \rho u_z\rangle }_{\varphi }^H-{\Omega }^{-H}{\langle \rho u_z\rangle }_{\varphi }^{-H}\right)+r\left({\langle \rho W_{z\varphi }\rangle }_{\varphi }^H-{\langle \rho W_{z\varphi }\rangle }_{\varphi }^{-H}\right)=0,\end{array}$$(C.5)

where we also define a vertical density-weighted average as $${\langle X\rangle }_H=\frac{1}{{\Sigma }_H\Delta r}{\displaystyle {\int }_{-H}^H{\displaystyle {\int }_{r-\Delta r/2}^{r+\Delta r/2}{\langle \rho \rangle }_{\varphi }}}X\text{d}r\text{d}z.$$(C.6)

The second term in the first line contains the radial transport due to the laminar flow and the stress tensors, while the second line corresponds to the same contributions across the disk surface.

Assuming a stationary disk, we ignored the time derivatives. Combining Eqs. (C.4) and (C.5) and assuming that the disk is symmetric with respect to the midplane (Ω^H = Ω^−H, we recall that $\Omega ={\langle u_{\varphi }\rangle }_{\varphi }/r$ by definition), we can derive $$\begin{array}{l}{\Sigma }_H{\langle u_r\rangle }_{\rho }\frac{\partial }{\partial r}\left(r^2{\Omega }^H\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r^2{\Sigma }_H{\langle W_{r\varphi }\rangle }_H\right)\\ +r\left({\langle \rho W_{z\varphi }\rangle }_{\varphi }^H-{\langle \rho W_{z\varphi }\rangle }_{\varphi }^{-H}\right)-\frac{1}{r}\frac{\partial }{\partial r}\left(r^3{\Sigma }_H{\langle \frac{\partial \Omega }{\partial z}{\displaystyle \int \rho }{u}_r\text{d}z\rangle }_{\rho }\right)=0.\end{array}$$(C.7)

The last two terms disappear when the integration is performed over the whole vertical extent with vanishing boundary conditions and the rotation profile is assumed vertically invariant.

The radial mass accretion rate within the disk $$\dot{M}(r)=-2\pi R{\Sigma }_H{\langle u_r\rangle }_{\rho }$$(C.8)

no longer has a simple expression. We limited our discussion to the vertical extent where the disk is practically Keplerian, and thus the last term in Eq. (C.7) can be neglected. We can then infer Eq. (14).

There are two contributions to the mass accretion due to angular momentum removal. The first is in the radial direction and the contribution is integrated across the disk thickness, for which we can define a mean α value such that $${\alpha }_r=\frac{{\langle W_{r\varphi }\rangle }_H}{{\langle c_{\text{s}}^2\rangle }_{\rho }}=\frac{{\widehat{T}}_{r\varphi }^H}{{\rho }_H{\langle c_{\text{s}}^2\rangle }_{\rho }},$$(C.9)

where ρ_H = Σ_H∕(2H) is the disk vertical mean density. The second contribution comes from the angular momentum that is lost due to vertical transport across both surfaces of the disk, and we define the corresponding α as $${\alpha }_z=\frac{(T_{z\varphi }^H-T_{z\varphi }^{-H})}{{\rho }_H{\langle c_{\text{s}}^2\rangle }_{\rho }}=\frac{({\langle \rho \rangle }_{\varphi }^H{W}_{z\varphi }^H-{\langle \rho \rangle }_{\varphi }^{-H}{W}_{z\varphi }^{-H})}{{\rho }_H{\langle c_{\text{s}}^2\rangle }_{\rho }}.$$(C.10)

The normalization is always performed with respect to the bulk-averaged pressure.

Appendix D Numerical discretizing effects on the evaluation of the stress tensor

The azimuthally averaged stress tensor of the velocity is defined locally in a narrow range of radius Δr around position r, such that rapid radial fluctuations are smoothed out, $$W_{p\varphi ,\text{loc}}=\frac{1}{ϵ}{\displaystyle {\int }_{r-ϵ/2}^{r+ϵ/2}{\langle \delta u_p\delta u_{\varphi }\rangle }_{\varphi }}\text{d}r,$$(D.1)

where ϵ is a minimum value such that the function is smooth. When numerically evaluating the stress tensor, this average is performed within a practically larger region. We discuss here the effects of discretized calculation on the stress tensor. When calculating the stress tensorin a larger finite region, there is an additional velocity difference with respect to the mean velocity in the whole region. We have thus $$\begin{array}{ll}{W}_{p\varphi ,\text{fin}}\hfill & =\frac{1}{\Delta r}{\displaystyle {\int }_{r-\Delta r/2}^{r+\Delta r/2}\langle (}\Delta u_p+\delta u_p)(\Delta u_{\varphi }+\delta u_{\varphi }){\rangle }_{\varphi }\text{d}r\hfill \\ \hfill & ={\overline{W}}_{p\varphi ,\text{loc}}+\frac{1}{\Delta r}{\displaystyle {\int }_{r-\Delta r/2}^{r+\Delta r/2}{\langle \Delta u_p\Delta u_{\varphi }\rangle }_{\varphi }}\text{d}r\hfill \\ \hfill & ={\overline{W}}_{p\varphi ,\text{loc}}+W_{p\varphi ,\text{cor}},\hfill \end{array}$$(D.2)

where Δ specifies the difference between the mean value at r and the mean value averaged in [−Δr∕2, Δr∕2], and $\overline{\cdot }$ stands for the radially averaged value. Below we perform a calculation with some simplifying assumptions to estimate the effect of averaging in a finite region on the radial component of the tensor.

Fig. D.1 Top: Reynolds α calculated atvarious bin sizes in r. The function fluctuates strongly at small bin sizes and becomes smooth for bin sizes ≥ 4d x. Bottom: gravitational α calculated atvarious bin sizes in r. The function fluctuates strongly between positive and negative values regardless of the increased bin size. It is impossible to measure the gravitational α in the simulations in this study, probably because of Cartesian discretization. The resolution dx = 0.12 AU.

First of all, we assumed Keplerian rotation such that the azimuthal average gives $${\langle u_{\varphi }\rangle }_{\varphi }(r)=r\Omega (r)=\sqrt{\frac{G{M}_{\ast }}{r}}.$$(D.3)

In the radial direction, the accretion velocity is given by $${\langle u_r\rangle }_{\varphi }(r)=-\frac{{(r^2\Sigma W_{r\varphi ,\text{loc}})}^{\prime }}{r\Sigma {(r^2\Omega )}^{\prime }}~-{\eta }_1\frac{W_{r\varphi ,\text{loc}}}{r\Omega }=-{\eta }_1{W}_{r\varphi ,\text{loc}}\sqrt{\frac{r}{G{M}_{\ast }}},$$(D.4)

where η₁ is a factor on the order of unity from the unknown profiles of Σ and W_rϕ,loc. By assuming the average inside Δr is the value at r, we can calculate $$\begin{array}{ll}\Delta u_r\Delta u_{\varphi }\hfill & =-{\eta }_1{W}_{r\varphi ,\text{loc}}\left(\sqrt{r+\Delta r/2}-\sqrt{r}\right)\left(\frac{1}{\sqrt{r+\Delta r/2}}-\frac{1}{\sqrt{r}}\right)\hfill \\ \hfill & =-{\eta }_1{W}_{r\varphi ,\text{loc}}\left(2-\frac{2r+\Delta r/2}{\sqrt{(r+\Delta r/2)r}}\right).\hfill \end{array}$$(D.5)

To calculate the average in the whole region, we assumed that W_rϕ,loc(r) and Σ(r) vary slowly and can be locally approximated as a constant. We also neglected the disk geometry and performed a linear integration along r. The correction term for averaged stress tensor calculated in a non-negligible finite region of [−Δr∕2, Δr∕2] is thus $$\begin{array}{ll}{W}_{r\varphi ,\text{cor}}\hfill & =\frac{1}{\Delta r}{\displaystyle {\int }_{r-\Delta r/2}^{r+\Delta r/2}{\langle \Delta u_r\Delta u_{\varphi }\rangle }_{\varphi }}\text{d}r\hfill \\ \hfill & =\frac{{\eta }_2}{\Delta r}{\displaystyle {\int }_{-\Delta r/2}^{\Delta r/2}\langle {\Delta }^{\prime }{u}_r{\Delta }^{\prime }{u}_{\varphi }\rangle }\text{d}{\Delta }^{\prime }r\hfill \\ \hfill & ={\eta }_1{\eta }_2{W}_{r\varphi ,\text{loc}}\left[\frac{1}{48}{\left(\frac{\Delta r}{r}\right)}^2-\frac{1}{128}{\left(\frac{\Delta r}{r}\right)}^3+\mathcal{O}{\left(\frac{\Delta r}{r}\right)}^4\right],\hfill \end{array}$$(D.6)

where η₂ is also a factoron the order of unity from the simplifications. The stress tensor that we measure from the simulations is therefore $$W_{r\varphi ,\text{fin}}={\overline{W}}_{r\varphi ,\text{loc}}+W_{r\varphi ,\text{cor}}\approx {\overline{W}}_{r\varphi ,\text{loc}}\left[1+\frac{{\eta }_1{\eta }_2}{48}{\left(\frac{\Delta r}{r}\right)}^2\right].$$(D.7)

As long as Δr∕r is small and W_rϕ,loc varies slowly, the correction term remains negligible and does not affect the estimation of the stress tensor. Nonetheless, the choice of Δr∕r still has to be reasonable such that the simplifying assumptions stay valid.

When the α values are measured from our simulations, the size of the binning has to be carefully chosen. As shown in Fig. D.1 (top panel), the value of the Reynolds α fluctuates strongly when the bin size in r is taken to be the smallest cell size. Nonetheless, by increasing the bin size, the measured α converges toward a smooth function. This function does not seem to be dependent of the bin size as long as the bin size ≥ 4dx, confirming the above discussions.

On the other hand, this does not seem to be true for the gravitational α. The measurement remains noisy regardless of the increased bin size. This might be a result of noise introduced by the Cartesian discretization for the disk system.

References

All Tables

All Figures

Fig. 1 Column density snapshots of the canonical simulation R_ℓ14. Left: edge-on. Right: face-on. The coordinates are in AU, and the origin is at the box center. The images are centered on the sink particle, and the white dashed circle indicates the numerical radius of sink accretion (4dx). After the formation of the sink particle, the originally roundish structure of the flattened first Larson-core is quickly accreted and a flat disk forms, with sightly smaller size. The protoplanetary disk does not stay at the center of the computation box and might experience a complex accretion history. The size of the disk grows slowly from its initial radius of about 20 AU and presents fluctuations. At about 120 kyr, the disk grows more rapidly and reaches almost 100 AU. Although only a few snapshots are presented here, the actual picture is highly variant in time.

In the text

	Fig. 2 Zoomed views of the disk at age (103.1 kyr in simulation time) from edge-on (left column) and face-on (right column). Top row: R_ℓ14. Bottom row:R_40ky_ℓ18. The column density is shown as in Fig. 1. The disk interior is clearly more resolved in R_40ky_ℓ18 with 0.06 AU resolution, particularly in the vertical direction.
In the text

Fig. 3 Upperpanel: source function as proposed by Ulrich (1976) and Hueso & Guillot (2005). The radius is normalized to the centrifugal radius rc, and the source term is normalized to the mean surface mass flux (${\dot{M}}_{\text{tot}}/(\pi r_{\text{c}}^2)$) such that its surface integral gives unity. Lower panel: cumulated mass accretion rate (normalized by the total mass accretion rate), Ṁ, integrated within radius r. The values of Ṁtot and rc are to be specified as functions of time according to the initial conditions of the prestellar core envelope.

In the text

Fig. 4 Scale heights from vertical hydrostatic equilibrium, hp (blue), from the first (h1, orange) and second (h2, green) moment of mass, plotted against the radius. Top two panels: R_ℓ14 and R_40ky_ℓ18 at 103 kyr. Bottom two panels: R_ℓ14 and R_80ky_ℓ18 at 138 kyr. The horizontal and the left vertical gray lines mark the numerical resolution (dx, outside the figure extent, thus not shown for the high-resolution runs). The values are averaged within rings of thickness dx and plotted against the mean radius of the rings. The second vertical gray line marks the radius of the numerical sink accretion scheme (4dx). The two vertical dotted lines mark the characteristic radii of the disk (Rkep and Rmag), as we defined in Sect. 5.1.2. Within Rkep the disk is close to vertical thermal equilibrium and the three definitions of scale height coincide very well. Beyond this radius, magnetic support becomes strong and deviation occurs. The thermal scale height, hp (r), is fitted to apower law, which is shown in the legend (in units of AU).

In the text

Fig. 5 Disk surface density profile Σcyl(r) measured inside a selected cylindrical region (blue). Top two panels: R_ℓ14 and R_40ky_ℓ18 at 103 kyr. Middle: radially accumulated mass of the disk of R_40ky_ℓ18 at 103 kyr. Bottom two panels: R_ℓ14 and R_80ky_ℓ18 at 138kyr. A three-segment power-law fit is overplotted (orange), and the power-law exponentsare shown in the legend. Two vertical dotted lines show the characteristic radii Rkep and Rmag that correspond to the transition of the power-law slope. The masses of the star M* and of the disk Md (inside Rkep plus betweenRkep and Rmag) are also displayed. Most of the disk mass is contained within Rkep.

In the text

Fig. 6 Radius evolution of the disk. The radii measured from simulations as described in Sect. 5.1.3 are plotted with solid (Rkep) and dashed (Rmag) lines for R_ℓ14 in blue. The disk initially has a Keplerian radius of about 20 AU for almost 50 kyr and a magnetized envelope that is roughly 1.5 times larger. At later times, the disk grows to nearly 50 AU in radius, and the magnetized region grows even more significantly to more than twice the size of the demagnetized Keplerian disk. Runs R_40ky_ℓ18 and R_80ky_ℓ18 are overplotted (with the same line styles), and the size of the disk is consistent with different resolutions.

In the text

Fig. 7 Mass evolution of the star-disk system. Run R_ℓ14 is plotted as a thick blue line for the sink mass, a thin line for the mass within Rkep, and as a dashed line for the mass within Rmag. A flattened core-disk structure forms at the beginning. After the sink forms, it quickly accretes most of the mass inside the disk, and the disk mass drops to about 1% of the stellar mass. The magnetized region around the Keplerian disk, though significant in size, contains much less mass. The high-resolution restarts are overplotted with the same line styles, and there is an overall good agreement.

In the text

Fig. 8 Mass accretion rate evolution. The canonical run at 1 AU resolution is plotted in blue. The mass accretion rate of the sink particle (thick line) is derived directly from its mass evolution. The accretion rate onto the disk is the integral of the measured source function (dashed line, as described in Sect. 5.5.3). The mass accretion rate is initially high an decreases slowly with time. The mass accretion onto the sink particle, though generally decreasing in time, shows some bursty behavior. This is due to the accumulation at the inner disk border shortly before reaching the accretion threshold. After 130 kyr, the density of the disk drops so much such that the threshold density of the sink accretion algorithm is no longer met and the sink hardly accretes. The high-resolution restarts are overplotted (solid line for sink accretion and stars for disk surface accretion). The mass flux onto the disk surface is unchanged in general, while the sink accretion is sensitive to the inner boundary condition of the disk (i.e., the choice of nacc, see Sect. 3.3).

In the text

Fig. 9 Radial temperature profiles of the disk. The same snapshots as in Fig. 4 are shown. The temperature is averaged (mass-weighted) vertically within H (orange) and3H (green). The snow-line (R150) and the CAI condensation line (R1500, if present) are marked with vertical dashed cyan lines. In the low-resolution run, the snow-line R150 is very close to the resolution limit and thus not well resolved. High-resolution restarts are indeed necessary to study the interior structure of the disk correctly. The CAI condensation line R1500 appears in R_40ky_ℓ18, while in R_80ky_ℓ18, this line hasmigrated inward due to the decreased overall density and temperature.

In the text

Fig. 10 Radial temperature profiles of the disk with and without stellar irradiation at time 120 kyr. Values of R_ℓ14_nofeed are plotted with lighter colors. The temperature is averaged (mass-weighted) vertically within H (orange) and3H (green). The snow line (R150) is marked with a vertical dashed cyan line. The disk temperature is significantly lower when no stellar irradiation is introduced, while the slope of radial profile does not seem to differ significantly.

In the text

Fig. 11 Evolution of the water snow line, R150. The snow line is initially located at about 10 AU and migrates inward in time as the disk evolves and loses its mass, finally stabilizing at ~5 AU. High-resolution runs are overplotted. R_80ky_ℓ18 has consistent temperature, while in R_40ky_ℓ18 the snow line migrates to a larger radius, possibly due to the mass accumulation in the disk and the increase in opacity. This reflects the importance of properly choosing nacc with increased resolution (see Sect. 3.3). The run with no stellar irradiation is also plotted for comparison, where the temperature is significantly lower.

In the text

Fig. 12 Velocity profiles of the disk at two resolutions: rotational velocity within H (orange) /3H (green) and infall velocity within H (red)/3H (violet). Both positive and negative values with |u| > 102 are plotted in logarithmic scale, and this is connected with a linear range in between. The Keplerian rotational velocity is plotted in blue for reference. The inner disk is very close to Keplerian rotation, and the deviation becomes remarkable in the outer magnetized part of the disk. The envelope is far from Keplerian rotation, while the infalling gas has roughly half the free-fall velocity. Below the disk accretion radius (vertical gray line), the velocity is not physically meaningful. The disk characteristic radii Rkep and Rmag (the dotted vertical lines) mark the transition between different behaviors.

In the text

Fig. 13 Plasma beta, the ratio between thermal and magnetic pressures, β = Ptherm∕Pmag, evaluated within one (orange) and three (green) times the disk scale height. The disk is highly demagnetized because ambipolar diffusion leaves the field lines outside during the collapse. The β value passes through unity around the junction between the disk and the envelope. The difference between the two curves comes from the lower magnetization level close to the disk midplane.

In the text

	Fig. 14 Disk α for R_40ky_ℓ18 at 103 kyr. Top panel: Reynolds and Maxwell stresses in r direction. Bottom panel: two components in z direction. Both positive and negative values with |α| > 10−4 are plotted in logarithmic scale, and this is connected by a linear range in between.
In the text

Fig. 15 α evolution of R_ℓ14 (increasing time from left to right). The four components αRey,r, αMax,r, αRey,z, αMax,z (from top to bottom). The r components are vertically averaged values between zero and a given altitude, and the z components are local values. The white curves trace the thermal scale height H and 3H. Positive values (red) can be roughly interpreted as a radial accretion, and negative values (blue) as radial expansion. The resolution is 1 AU. As the disk grows in size, the interior of the disk is described with more cells and there is a clearer pattern of α. The early disk has relatively strong αRey compared to its magnetic counter part. It is probably under constant perturbation by the strong infall from the envelope. The first two chosen snapshots are not necessarily representative of this epoch because the behavior of αRey is very irregular in time, leading to a disk that fluctuates in size around 10 AU. The infall from the envelope weakens later, and the disk becomes more regular, with its inner parts readily accreting. This accretion is mostly caused by vertical magnetic transport of angular momentum within the envelope.

In the text

Fig. 16 Measured α values of R_ℓ14 (left) and R_40ky_ℓ18 (right) at 103 kyr. The disk radius is about 15 AU. The bins for averaging are chosen to be larger than dx in the high-resolution run in order to remove numerical noise. These figures clearly show that the technique of restart using higher resolutions is needed to study the interior properties of the disk. The z components are more important because they dominate in governing the mass accretion inside the disk. The middle layers of the disk are more irregular with some partially expanding parts, while most of the accretion occurs through the upper layer slightly above 3H and reaches a small radius before penetrating the surface of the disk.

In the text

Fig. 17 Measured α values of R_ℓ14 (left) and R_80ky_ℓ18 (right) at 138 kyr. Disk radius (Rkep) grows to 50 AU at this time. The values are evaluated on a grid of 1 AU resolution for both runs. The global patterns are compatible for the two runs with 16 times difference in resolution. The middle layers of the disk are radially structured and partially expending at times, while the accretion occurs through a channel slightly above 3H and brings a significant amount of mass from the infalling envelope to a small radius without passing through the disk.

In the text

	Fig. 18 Toomre Q ≈ Ωcs∕(πGΣ) parameter evaluated within one (orange) and three (green) times the disk scale heights. The same snapshots are taken as in Fig. 4. The disk is very stable against self-gravitating fragmentation. The difference between the two curves mainly arises because the total mass is not taken into account for the value within H and 3H.
In the text

Fig. 19 Poloidal mass flux (${\langle \rho u_p\rangle }_{\varphi }$) and specific angular momentum ($r{\langle u_{\varphi }\rangle }_{\varphi }$) around the disk region (color streamlines in the gray image, in SI units), for R_ℓ14 (top) and R_40ky_ℓ18 (bottom)at simulation time 103 kyr. All quantities are azimuthally averaged and averaged between the upper and lower planes. White curves outline the disk scale height H and 3H at r < Rkep. An outflow is launched at high altitudes at z ≈ 20 AU in R_ℓ14. In the high-resolution run, the flow varies slightly more, but it is qualitatively similar. The outflow is launched closer to the disk (z ≈ 10 AU) and the disk size is slightly smaller. The outflow has relatively low density and manifestsitself as a cavity. Significant infall from the envelope is channeled through a zoneof low angular momentum at several scale heights (darker gray region near z ≈ r) of the disk. The horizontal infall from the envelope reaches the outer edge of the disk (shown in Fig. 2 at r ≈ 15 AU), and thedensity contrast generates a shock that pushes the gas to higher altitude. This creates a counterclockwise poloidal circulation near the disk edge between 10 and 15 AU in R_40ky_ℓ18 because the outer disk decretes at this instant, while it is not seen in the upper panel.

In the text

	Fig. 20 Zones showing the main properties of the disk-envelope system in Sect. 5.5.1 and 5.5.2.
In the text

Fig. 21 Poloidal mass flux (${\langle \rho u_p\rangle }_{\varphi }$) and specific angular momentum ($r{\langle u_{\varphi }\rangle }_{\varphi }$) around the disk region (color streamlines in the gray image, in SI units) for R_ℓ14 (top) and R_80ky_ℓ18 (bottom) at simulation time 138 kyr. All quantities are azimuthally averaged and averaged between the upper and lower planes. The white curves outline the disk scale height H and 3H at r < Rkep. An outflow is launched at z ≈ 20 AU. The disk has grown in size, while the flow pattern stays very similar to that in Fig. 19. At this moment, the magnetized zone is almost twice the size of the Keplerian disk (Rmag ≈ 2Rkep), and the counterclockwise poloidal circulation around Rkep is evident inboth panels and occupies almost the entire magnetized zone and the outer Keplerian disk.

In the text

Fig. 22 Angular momentum flux and angular momentum ($r{\langle \rho u_{\varphi }\rangle }_{\varphi }$) around the disk region (color streamlines in the gray image, in SI units), for R_ℓ14 (left) and R_40ky_ℓ18 (right) at 103 kyr. All quantities are azimuthally averaged and averaged between the upper and lower planes. From top to bottom: Laminar transport ($r{\langle u_{\varphi }\rangle }_{\varphi }{\langle \rho u_p\rangle }_{\varphi }$), turbulent transport ($r{\langle \rho \delta u_{\varphi }\delta u_p\rangle }_{\varphi }$), magnetic transport ($-r{\langle B_{\varphi }{B}_p\rangle }_{\varphi }/(4\pi )$), and their sum. White curves outline the disk scale height H. The laminar flow convects mass with positive angular momentum, while the tensor terms (Reynolds and Maxwell) transport the angular momentum outward. In particular, the magnetic tensor dominates in the envelope, and this allows the mass to fall onto the disk through the loss of angular momentum. The general pattern persists with improved resolution. White patches in the gray image are due to locally negative angular momentum. This is an artifact of the averaging because the selected axis is not exactly aligned with the asymmetric outflow.

In the text

Fig. 23 Angular momentum flux and angular momentum ($r{\langle \rho u_{\varphi }\rangle }_{\varphi }$) in the disk surrounding (color streamlines in the gray image, in SI units) for R_ℓ14 (left) and R_80ky_ℓ18 (right) at 138 kyr. The disk scale height H is shown in white. All quantities are azimuthally averaged and averaged between the upper and lower planes. From top to bottom: Laminar transport ($r{\langle u_{\varphi }\rangle }_{\varphi }{\langle \rho u_p\rangle }_{\varphi }$), turbulent transport ($r{\langle \rho \delta u_{\varphi }\delta u_p\rangle }_{\varphi }$), magnetic transport ($-r{\langle B_{\varphi }{B}_p\rangle }_{\varphi }/(4\pi )$), and their sum.

In the text

Fig. 24 Source function S(r) measured at H (orange) and3H (green) in R_ℓ14 and R_40ky_ℓ18 at 103 ky, R_ℓ14 and R_80ky_ℓ18 at 138 ky (from top to bottom). The flux measured at the horizontal plane z = 4d x is shown in gray. The curves in lighter colors are the analytical functions from pure hydrodynamic predictions that have the same amount of total flux (see Fig. 3). The source function measured in the simulations is much steeper than that of the classical model, and strong infall occurs in the central part of the disk. At a given radius, the flux is generally greater at lower altitudes, indicating that there is a strong radial flux in the interior of the disk. At small radii (below the sink accretion radius), the measured flux at H∕3H is no longer realistic, and we used the measurement at z = 4dx to calculate the total flux.

In the text

Fig. 25 Massaccretion rate measured across the disk surface (Ṁs(< r)), measured radially inside the disk (Ṁr(r)), and inferred from the measured α (Ṁα (r)) for R_ℓ14 and R_40ky_ℓ18 at 103 ky, R_ℓ14 and R_80ky_ℓ18 at 138 ky (from top to bottom). Both positive and negative values with |Ṁ | > 10−7 are plotted in logarithmic scale, and this is connected with a linear range in between. Most of the surface flux arrives where the slope of Ṁs(< r) is the steepest. Below this radius, the disk readily accretes in the radial direction only at a very small radius (<1 AU), and the outer part can have complex motions. The radial mass flux derived from α is not exactly the same as Ṁr(r), although the amplitudes are comparable.

In the text

	Fig. B.1 Temperature plotted against altitude, z, at various distance to the central star, r, from R_ℓ14 (the disk aspect ratio is ~0.1 (Sect. 5.1.1). Every dot corresponds to one cell in the simulation, without any averaging. The discrete distribution comes from the slight difference between the grid direction and the selected disk axis.
In the text

Fig. D.1 Top: Reynolds α calculated atvarious bin sizes in r. The function fluctuates strongly at small bin sizes and becomes smooth for bin sizes ≥ 4d x. Bottom: gravitational α calculated atvarious bin sizes in r. The function fluctuates strongly between positive and negative values regardless of the increased bin size. It is impossible to measure the gravitational α in the simulations in this study, probably because of Cartesian discretization. The resolution dx = 0.12 AU.

In the text