Early evolution of disk, outflow, and magnetic field of young stellar objects: Impact of dust model

The formation and early evolution of low mass young stellar objects (YSOs) are investigated using three-dimensional non-ideal magneto-hydrodynamics simulations. We investigate the evolution of YSOs up to ~ 10^4 yr after protostar formation, at which protostellar mass reaches ~ 0.1 M_\odot . We particularly focus on the impact of the dust model on the evolution. We found that a circumstellar disk is formed in all simulations regardless of the dust model. Disk size is approximately 10 AU at the protostar formation epoch, and it increases to several tens of AU at ~ 10^4 yr after protostar formation. Disk mass is comparable to central protostellar mass and gravitational instability develops. In the simulations with small dust size, the warp of the pseudodisk develops ~ 10^4 yr after protostar formation. The warp strengthens magnetic braking in the disk and decreases disk size. Ion-neutral drift can occur in the infalling envelope under the conditions that the typical dust size is a \gtrsim 0.2{\mu}m and the protostar (plus disk) mass is M \gtrsim 0.1 M_\odot. The outflow activity is anti-correlated to the dust size and the strong outflow appears with small dust grains.


INTRODUCTION
Molecular cloud cores, which are birth places for protostars and protoplanetary disks, are strongly magnetized. Measurement of the Zeeman effect has shown that the massto-flux ratio normalized by its critical value is of order unity (Troland & Crutcher 2008;Crutcher 2012). The strong magnetic field of cloud cores is also supported by threedimensional simulations of molecular cloud formation (Inoue & Inutsuka 2012). This suggests that although the magnetic field is not strong enough to support the core against gravitational collapse, it should play a crucial role during a gravitational collapse of the core. For example, angular momentum removal from the central region by the magnetic field (so called magnetic braking) almost completely suppresses disk formation if the neutral gas and magnetic field are well coupled (e.g., Mellon & Li 2008;Tsukamoto et al. 2015b).
Another important feature of cloud cores is their low ionization degree (e.g., Umebayashi & Nakano 1990;Nakano et al. 2002). The ionization degree of cloud cores (∼ 10 −20 g cm −3 ) is typically ∼ 10 −8 , and it decreases as den-sity increases during the gravitational contraction phase. In such a weakly ionized magnetized cloud, non-ideal effects, ohmic diffusion, the Hall effect, and ambipolar diffusion play key roles.
It has been shown that ohmic diffusion largely affects the formation of young stellar objects. Ohmic diffusion decouples the gas from the magnetic field at ρ > 10 −11 g cm −3 (Umebayashi & Nakano 1990;Nishi et al. 1991;Nakano et al. 2002) and gas can accrete to the central protostar leaving the magnetic flux at this density. Furthermore,  and Tsukamoto et al. (2015b) showed that disk formation at the protostar formation epoch with size of r ∼ 1 to 10 AU is enabled by ohmic diffusion due to the decoupling. However, note that the magnetic flux accumulates around the central dense region of ρ 10 −11 g cm −3 in the protostellar accretion phase (the evolution phase after protostar formation), and a huge amount of magnetic flux abides in the central compact region when we neglect other non-ideal effects; this is because ohmic resistivity is an increasing function of density and does not depend on the magnetic field strength. Thus, neglecting other non-ideal ef-fects may affect evolution after protostar formation, during which a large amount of magnetic flux is supplied towards the central region.
Ambipolar diffusion, on the other hand, plays a role not only in the high-density region but also in the low-density region of ρ 10 −11 g cm −3 , and it becomes strong as magnetic flux accumulates. Therefore it must play a central role in magnetic field evolution in the protostellar evolution phase. Previous studies have shown that ambipolar diffusion further weakens the coupling between the magnetic field and the gas in a newly born disk (Tsukamoto et al. 2015b). Furthermore, magnetic flux can drift outwardly relative to the neutral motion in the envelope by ambipolar diffusion (e.g., Li et al. 2011;Tsukamoto et al. 2017a;Zhao et al. 2018b). This outward magnetic field drift is a promising mechanism to remove the magnetic flux from the central region in the protostellar accretion phase. Furthermore, magnetic field drift possibly accompanies ion drift in the low-density region, or more precisely, in the region where Hall parameter is β Hall 1. This may provide a unique observational opportunity to quantify the magnetic field in the inner envelope of YSOs (see for recent attempt by Yen et al. 2018).
Because non-ideal effects (or finite conductivity) arise due to the low ionization degree of the cloud core, they inevitably depend on the microscopic properties of the gas. Previous studies have shown that cosmic-ray ionization and dust size distribution are the keys to determining the resistivity of non-ideal effects (Nishi et al. 1991;Zhao et al. 2018a; Koga et al. 2019). They are the source and sink of charged particles, respectively.
There are observations suggesting that dust growth may proceed even at cloud core scale. For example, mid-infrared emission from the cloud core is interpreted as scattered light of µm sized dust grains Pagani et al. 2010). In theoretical estimates of dust growth, dust size of 1µm is possibly realized within the free-fall time in the dense inner part of the cloud core (Ormel et al. 2009;Hirashita & Li 2013). Thus, the impact of dust size on the non-ideal effect and cloud core evolution should be investigated.
Recently, Zhao et al. (2016Zhao et al. ( , 2018b suggested that dust growth and removal of the small dust grains from MRN dust size distribution ("truncated MRN") changes magnetic resistivity significantly and enables disk formation. On the other hand, several studies have reported disk formation with small sized dust grains Wurster et al. 2016;Tomida et al. 2015;Tsukamoto et al. 2017a;Wurster et al. 2018;Wurster & Bate 2019). Therefore, there are inconsistencies among previous studies. Zhao et al. (2018b) employed several simplifications for numerical treatments to avoid small time stepping such as setting upper limit on ambipolar resistivity and Alfvèn velocity, neglecting Ohmic diffusion, and imposing relatively large inner boundary of 2 AU. Note that large sink radius of 1 AU tends to suppress disk formation . We speculate that these treatments may have some impacts on disk formation and early evolution.
Thus, we believe that an additional study investigating the impact of dust models is required. In this study, we investigated the early disk evolution phase up to ∼ 10 4 yr after protostar formation particularly focusing on the impact of dust size difference.

Numerical Method
In our numerical simulations, non-ideal magnetohydrodynamics (MHD) equations were solved.
where ρ is gas density, P is gas pressure, B is the magnetic field,B is defined asB ≡ B/|B|. ηO and ηA are the resistivity for ohmic and ambipolar diffusion, respectively, and Φ is the gravitational potential. G is a gravitational constant. We adopted a barotropic equation of state (EOS) in which gas pressure only depends on density. cs,iso = 190 m s −1 is isothermal sound velocity at 10 K. We used a critical density of ρcrit = 4 × 10 −14 g cm −3 above which gas behaves adiabatically. In this study, we ignored the Hall effect. We used the smoothed particle magneto-hydrodynamics (SPMHD) method to solve the equations (Iwasaki & Inutsuka 2011. Our numerical code was parallelized with Message Passing Interface (MPI). We treated ohmic and ambipolar diffusion according to the prescription described in Tsukamoto et al. (2013a). Both the diffusion processes were accelerated by super-time stepping (STS) (Alexiades et al. 1996).
To calculate the time evolution after protostar formation, we employed the sink particle technique (Bate et al. 1995). The sink particle was dynamically introduced when the density exceeds ρ sink = 4 × 10 −13 g cm −3 . In our simulations, one sink particle was permitted. The sink particle absorbs SPH particles with ρ > ρ sink within r < r sink = 1 AU, and the mass and linear momentum of SPH particles are added to those of sink particle. The sink particle interacts with SPH particles via gravity. The system was integrated up to ∼ 10 4 yr after protostar formation, at which the mass of the sink particle (or protostar) reached ∼ 0.1 M .

Resistivity model
We used the tabulated resistivity calculated by the methods described in Susa et al. (2015). We considered the ion species of H + 3 , H + 2 , H + 3 , HCO + , Mg + He + , C + , O + , O + 2 , H3O + , OH + , H2O + and the neutral species of H, H2, He, CO, O2, Mg, O, C, HCO, H2O, OH. We also considered the neutral and singly charged dust grains, g 0 , g − , g + . We took into account cosmic-ray ionization, gas-phase and dust-surface recombination, and ion-neutral reactions. We also considered the indirect ionization by high-energy photons emitted by direct cosmic-ray ionization (described as CRPHOT in the UMIST database).
The initial abundance and reaction rates were taken from the UMIST2012 database (McElroy et al. 2013). The grain-ion and grain-grain collision rates were calculated using the equations of Draine & Sutin (1987). The chemical reaction network was solved using the CVODE package (Hindmarsh et al. 2005). We assumed that the system was in chemical equilibrium, which is valid in the nearby star forming cloud as discussed in Marchand et al. (2016). We calculated resistivity using the abundances of charged species in the equilibrium state. The momentum transfer rate between neutral and charged species was calculated using the equations described in Pinto & Galli (2008). The temperature for the chemical network was assumed to be T = 10(1 + γT (ρ/ρc) (γ T −1) ) K, where γT = 7/5. The dust internal density is fixed to be ρ d = 2 g cm −3 . The cosmic ray ionization is fixed to be ξCR = 10 −17 s −1 .
For the dust models, we considered single sized dust models and models with a size distribution of n(a) ∝ a −3.5 . For the single sized dust models, we considered three dust size of of a = 0.035µm, 0.1µm, and 0.3µm. For the models with size distributions, we considered models with minimum and maximum dust sizes of amin = 0.005µm and amax = 0.25µm (MRN distribution) and amin = 0.1µm and amax = 0.25µm (truncated MRN distribution which mimics the size distribution of Zhao et al. 2018b). We assumed a fixed dustto-gas mass ratio of 0.01.

Initial conditions
We adopted the density-enhanced Bonnor-Ebert sphere surrounded by medium with a steep density profile of ρ ∝ r −4 as the initial density profile, ρ(r) = ρ0ξBE(r/a) for r < Rc ρ0ξBE(Rc/a)( r Rc ) −4 for Rc < r < 10Rc.
and a = cs,iso f 4πGρ0 where ξBE is non-dimensional density profile of the critical Bonnor-Ebert sphere, f is a numerical factor related to the strength of gravity, and Rc = 6.45a is the radius of the cloud core. f = 1 corresponds to the critical Bonnor-Ebert sphere, and the core with f > 1 is gravitationally unstable. A Bonnor-Ebert sphere is determined by specifying central density ρ0, the ratio of the central density to density at Rc ρ0/ρ(Rc), and f . In this study, we adopted the values of ρ0 = 7.3×10 −18 g cm −3 , ρ0/ρ(Rc) = 14, and f = 2.1. Then, the radius of the core is Rc = 4.8×10 3 AU, and the enclosed mass within Rc is Mc = 1 M . The α therm (≡ E therm /Egrav) is equal to 0.4, where E therm and Egrav are the thermal and gravitational energy of the central core (without surrounding medium), respectively. The steep envelope was adopted to put the outer boundary far away from the central cloud core. With the steep profile, the total mass of the entire domain remains ∼ 2Mc. For rotation of the cloud core, we adopted an angular velocity profile of Ω(d) = x 2 + y 2 and Ω0 = 2.3 × 10 −13 s −1 . Ω(d) is almost constant for d < 1.5Rc and rapidly decreases for d > 1.5Rc. The ratio of the rotational to gravitational en-ergy βrot within the core is βrot (≡ Erot/Egrav) = 0.03, where Erot is the rotational energy of the core.
We constructed a magnetic field profile that is constant, has only the z component at the center, and asymptotically obeys B ∝ r −2 as r → ∞ (see Appendix A for details of the magnetic field configuration). The merit of this profile is the avoidance of a low β region in the surrounding medium, which appears when adopting constant magnetic field strength and radially decreasing density. We assumed the characteristic length scale of the field configuration to be R0 = Rc. The central magnetic field strength and plasma β were B0 = 62µG and β = 1.6 × 10 1 , respectively. With our magnetic field profile, the mass-toflux ratio of the core µ relative to the critical value was µ/µcrit = (Mc/Φmag) = 8, where Φmag is the magnetic flux of the core and µcrit = (0.53/3π)(5/G) 1/2 . The mass-to-flux ratio is large compared to the observed value because the magnetic field becomes weak in the outer region with our magnetic field configuration. However, note that if we assume a constant magnetic field profile with central magnetic field strength which is widely used in previous studies, the mass-to-flux ratio is µconst/µcrit = 4. This value would be more suitable to compare our magnetic field strength with those of previous studies because we investigated time evolution of the gas in the central region of the core. We resolved 1 M with 3 × 10 6 SPH particles. Thus, each particle had a mass of m = 3.3 × 10 −7 M .
The model names and corresponding dust sizes are summarized in Table 1. Table 1 also summarizes simulation results.

Impact of dust size on resistivity
First, we investigate how resistivity depends on the dust model. Figure 1 shows ηO and ηA with different dust models. The left panel of figure 1 shows ηO and ηA as a function of density. This indicates that an increase in dust size has two contradicting impacts on ηA. With larger dust size, ηA is large in the low-density region (ρ 10 −13 g cm −3 ) and small in the high-density region (ρ 10 −13 g cm −3 ). Thus, dust growth has a positive impact on the decoupling between magnetic field and gas in the low-density region but has a negative impact in the high-density region.
The small ηA with the small dust grains in the lowdensity region is caused by an increase of ion abundance in the gas phase. With small dust grains, the electrons in the gas phase are more efficiently absorbed by dust grains and ions lose their counterparts. As a result, recombination rate in the gas phase decreases. This causes an increase in ion abundance and decrease of ηA in the low-density region with small dust grains. Note that, as dust size increases, ηA becomes similar to the analytic formula by Shu (1983), ηA = B 2 /(γCρ 3/2 ) (dotted line) in the low-density region where γ = 3.5 × 10 13 cm 3 (g s) −1 and C = 3 × 10 −16 (g cm −3 ) 1/2 . This is because the abundance and total surface area of dust decrease as the dust size increases, and the system becomes similar to that of dust-free case. The mutually contradicting impact of dust size on ηA in the high and low-density regions may introduce diversity to the evolution of disk, outflow, and magnetic flux.  [AU] are the time at protostar formation, time at the end of the simulation, centrifugal radius at t end , outflow size at t end , and maximum radius at which magnetic field drift velocity is larger than 0.19 km s −1 at t end , respectively. On the other hand, ηO monotonically decreases as dust size increases. This is particularly clear from comparing the monosize dust model. This shows that ohmic diffusion becomes less important as dust grows in the molecular cloud core or in the circumstellar disk.
The right panel of figure 1 shows ηA as a function of the magnetic field at ρ = 4 × 10 −15 g cm −3 . This shows that the dependence of ηA on magnetic field strength is not simple particularly for small dust size. ηA of a=0.035µm is almost independent of the magnetic field in B < 10 −2 G (and hence behaves like "ohmic diffusion"). As the magnetic field increases, on the other hand, ηA obeys ηA ∝ B 2 . The change in dependence on magnetic field is related to the weak coupling between charged particles and the magnetic field at this density and magnetic field strength (for more detail, see Nakano et al. 2002). As dust size increases, the plateau slides to higher magnetic field strength (see blue and black lines) and becomes narrow. The right panel clearly indicates that, although the relation of ηA ∝ B 2 is often assumed for ambipolar diffusion, this is not always valid.

Time evolution of fiducial models
In this section, we describe the time evolution of two fiducial models: one with small dust grains and one with large dust grains. Important features, which will be discussed in subsequent sections, will be highlighted here.

Time evolution of a model with small dust grains
First, we investigate model MRN as a fiducial model with small dust grains. Figures 2 and 3 show the density evolution in the 500 AU scale box and in the 2000 AU box of model MRN, respectively. In this model, the protostar is formed at t = 4.1 × 10 4 yr. The top left panel of figure 2 shows that a weak outflow with v < 1 km s −1 is formed. This outflow expands to 150 AU in size at t = 4.5 × 10 4 yr (top middle) but almost stalls (top left). Then a stronger and more collimated outflow with v > 1 km s −1 is launched (bottom left). At t = 5.1 × 10 4 yr (which corresponds to ∼ 10 4 yr after protostar formation), fast outflow dominates in the envelope (bottom middle).
The slow outflow is driven by the first core and the fast outflow is driven by the circumstellar disk. In our simulations, the stellar outflow is not resolved due to the sink particle with radius of 1 AU. The outflow velocity is roughly determined by rotation velocity (i.e., Keplerian velocity; vK ∼ 1 km s −1 (100/AU) −1/2 (Mstar/0.1M ) 1/2 ) at the launching point. This suggests that the fast outflow is launched at several tens of AU corresponding to the disk radius.
The bottom middle panel of figure 3 shows that the outflow head reaches ∼ 1000 AU at this epoch. The property of outflow will be investigated in more detail in §3.8. We find that the warp of the pseudo-disk develops at r ∼ 100 AU, approximately 10 4 yr after protostar formation (bottom right panel). This pseudo-disk warp strengthens the magnetic field in the disk and negatively impacts disk growth. We will revisit pseudo-disk warp in more detail in §3.4. Figure 4 shows the plasma β map on the x-z plane in the 500 AU scale box. The plasma β has a large dynamic range from β > 10 3 in the disk to β < 10 −2 in the upper envelope. At the protostar formation epoch (top left panel), the low β region is localized in r < 100 AU, and it expands as time proceeds.
The white arrows around the midplane (around the x axis) indicate that the magnetic field is highly pinched towards the center, and the so-called "hour-glass" magnetic field configuration is realized. The white arrows in the bottom right panel show that the "neck" of the hour-glass magnetic field configuration shifts towards z < 0 and contracts as the warp develops. The warp of the pseudo-disk is more clearly seen in this β map. This contraction enhances the magnetic field in the disk. As a result, the plasma β of the disk decreases from β 10 3 to β ∼ 10 2 once the warped pseudo-disk develops (from bottom middle to bottom right panels). Figure 5 shows the density map on the x-y plane at the same epochs of figure 2 in the 250 AU scale box. The central high-density region (ρ 10 −14 g cm −3 ) is the circumstellar disk. In our simulations, the circumstellar disk forms immediately after protostar formation (or sink particle creation) and survives for t ∼ 10 4 yr thereafter. Thus, our results are consistent with the disk formation scenario in  and Inutsuka (2012) where the first core directly becomes the circumstellar disk. The disk radius gradually grows from ∼ 10 AU (top left) to ∼ 50 AU size (bottom middle). As the disk size increases, spiral arms develop due to gravitational instability. We observed that the spiral arms are repeatedly formed in the disk. The emergence of the spiral arms and outflow activity are correlated, and strong outflow is launched when the spiral arms are prominent. As the warp of the pseudo-disk develops, the disk begins to shrink (from bottom middle to right panels). More rigorous analysis of disk size evolution is presented in §3.5.

Time evolution of a model with large dust grains
Next, we will investigate the time evolution of a model with large dust grains, model a03µm. Figures 6 and 7 show the density evolution on the x-z plane in the 500 AU and 2000 AU scale boxes of model a03µm, respectively. The clear difference between this model and model MRN at the protostar formation epoch is the absence of outflow from first core (top left panels). This is due to stronger ambipolar diffusion in the low-density region (see figure 1). However, as time proceeds, the outflow with v > 1 km s −1 eventually forms at ∼ 7 × 10 3 yr after protostar formation (bottom left panels). The bottom left panel of figure 7 shows that the outflow is slightly tilted from the z direction. This is due to the non-axisymmetric spiral arms in the disk. (see figure 9; non-axisymmetric m=1 mode frequently develops in this model). In this model, the outflow is more collimated and its cylindrical radius ( x 2 + y 2 ) is smaller than that of model MRN at t = 5.3 × 10 4 yr. As shown §3.8, the outflow angular momentum of this model is much smaller than that of model MRN and the smaller cylindrical radius is one   reason for this difference. Another important difference between these two models is the absence of pseudo-disk warp. Even though we calculated the system evolution ∼ 10 4 yr after protostar formation, the pseudo-disk warp (and any symptom of it) does not develop in this model. Note that, in the case of model MRN, a weak warp is already formed at t = 5.1 × 10 4 yr (bottom middle panel of figure 2 and 4). Due to the absence of pseudo-disk warp, the large disk is maintained (see figure 9). Figure 8 shows the plasma β map of model a03µm. The top left panel shows that a bipolar low β structure already forms at the protostar formation epoch, but it does not drive the outflow until the magnetic field is sufficiently amplified by disk rotation. An interesting difference between model MRN and model a03µm is the thickness of the current sheet at the midplane. The bottom left panel shows that the contours of plasma β at r 50 AU are sparse around the midplane, indicating that the magnetic field slowly changes towards the vertical direction. Furthermore, the white arrows indicate that the magnetic field is weakly pinched towards the center. On the other hand, the bottom left panel of figure 4 show that the contours of plasma β in model MRN are dense around the midplane, and the white arrows show that the magnetic field is strongly pinched towards the center. This difference comes from the strength of ambipolar diffusion and significant magnetic field drift in the pseudodisk of model a03µm. Figure 9 shows the density evolution on the x-y plane. In this model, the circumstellar disk also forms immediately after protostar formation and survives for ∼ 10 4 yr after protostar formation. At its formation epoch, the disk size is similar to that of model MRN. As it grows, the radius becomes larger than that of model MRN. This is particularly clear in later epochs (bottom panels). The spiral arms are also more prominent in this model. We can clearly see that the disk monotonically grows in this model.

Magnetic field drift induced by ambipolar diffusion
One of the most important phenomena caused by ambipolar diffusion is magnetic field drift in the envelope. Magnetic field drift determines the magnetic field strength of a newly born circumstellar disk. It has also been suggested as a mechanism for magnetic flux redistribution in the envelope (e.g., Li 1998). Furthermore, ion-neutral drift, which accompanies magnetic field drift in the low-density region, may provide possible direct observational evidence of a relatively strong magnetic field in the envelope (see recent attempt of Yen et al. 2018). Figure 10 shows the azimuthally averaged radial velocities on the x-y plane at t = 5.3 × 10 4 yr (∼ 1.1 × 10 4 yr after protostar formation). The solid, dashed, and dotted lines show the gas radial velocity vr, radial drift velocity of magnetic field v drfit,r , and radial velocity of magnetic field v drfit,r + vr, respectively. The drift velocity is calculated as Among our simulations, model a03µm and  Early evolution of disk, outflow, and magnetic field of young stellar objects: Impact of dust model 9  model trMRN show significant radial drift in relatively extended region with size of r > 100 AU. In figure 10, we also plot v dirft,r of model MRN as an example of the model with small drift velocity.
In model a03µm, the drift velocity is positive and reaches ∼ 1 km s −1 at r ∼ 100 AU (dashed line). As a result, total radial velocity of the magnetic field (vr + v drift,r ) is ∼ −0.5 km s −1 (dotted line) and much slower than the gas infall velocity (solid line). We define r drift as the maximum value of radius at which the radial magnetic field drift velocity is larger than 0.19 km s −1 on x-y plane. Here we use velocity threshold of 0.19 km s −1 which corresponds to the sound velocity of T = 10 K and approximates the sound velocity in envelope. In model a03µm, r drift is ∼ 700 AU at t = 5.3 × 10 4 yr.
The radial drift velocity becomes small as the typical dust size decreases. In model trMRN, the drift velocity is typically ∼ 0.4 km s −1 . On the other hand, r drift also extends to ∼ 700 AU in this model. In model MRN, notable outward radial drift can not be observed, and vr +v drift,r (dotted line) is almost identical to the gas infall velocity (vr: solid line). This clearly shows that outward radial drift of the magnetic field occurs only with relatively large dust grains (or absence of small dust grains) in early evolution phase of young stellar objects. This is mainly because the difference of ηA in the low density region of ρ 10 −14 g cm −3 .
One of the most important aspects of magnetic field drift by ambipolar diffusion is that magnetic field drift accompanies ion-neutral drift. The ion-neutral drift is possibly observable as a velocity difference between for example CO and HCO + . Nakano et al. (2002) showed that the ratio of drift velocity of ion (v drift,ion ) and magnetic field (v drift ) can be calculated as where ξ is a constant of order unity, and β Hall is the Hall parameter, where qi and mi are the charge and mass of the ion, respectively. mn is the mean neutral mass. σv i is the rate coefficient for collisional momentum transfer between ion and neutral. For the ion species with β Hall 1, the ratio of velocity becomes v drift,ion /v drift ∼ 1 and the ion moves with the magnetic field, giving ion-neutral relative velocities. Figure 11 shows the azimuthally averaged Hall parameter on the x-y plane of model a03µm and of model trMRN at t = 5.3 × 10 4 yr. Here, we assume that HCO + is a major ion species in the envelope, and we use its mass and charge to calculate the Hall parameter. Furthermore, in this figure, we assume σv i is constant and σv i = 1.8 × 10 −9 cm 3 s −1 , which is enough for our purpose here. Note that the Hall parameter of other ion species is not significantly different because the mass and charge of the ion species in the envelope do not differ significantly. Figure 11 shows that Hall parameter is β Hall > 1 in the region of r 100 AU, and the ions are expected to move with the magnetic field in the envelope. Therefore, we conclude that the ion is well coupled with the magnetic field in the flattened envelope  and the magnetic-field drift velocity in figure 10 can also be regarded as the ion drift velocity in r 100 AU.
The size of the region in which ion-neutral drift occurs is important quantity to observe ion-neutral drift. Figure 12 shows the time evolution of r drift as a function of Mstar + M disk for model a03µm and for model trMRN. Here Mstar is the mass of the sink particle and M disk is the disk mass which is defined as enclosed mass of the region with ρ > 10 −14 g cm −3 (see §3.6). The figure shows that r drift monotonically increases as the mass (and hence magnetic flux) accumulates to the central region and reaches ∼ 100 AU at Mstar + M disk ∼ 0.1 M . This indicates that Mstar +M disk 0.1 M is required for magnetic field drift to occur in r > 100 AU and indicates that a sufficient amount of mass (and hence magnetic flux) should have accreted to the central region for magnetic field drift to occur in a relatively extended region. We discuss the interpretation of our results and the relation to the observations in §4.2.

Warp of pseudo-disk
As shown in the figures 2 and 4, an interesting structure develops in the pseudo-disk i.e., the warp of the pseudo-disk. Among our simulations, the warp appears in model MRN, model a01µm and model a0035µm at ∼ 10 4 yr after protostar formation. Figure 13 shows time evolution of density for model a01µm from t = 4.85 × 10 4 yr (just before the warp develops) to t = 5.1 × 10 4 yr and that how warp develops. At t ∼ 4.85 × 10 4 yr (top left panel), the density structure is approximately symmetric with respect to the x axis. Then, the magnetic field around the disk is perturbed by the spiral arm, and the neck of the hour-glass magnetic field is shifted to z < 0 direction (top right panel). The gas accretion flow is also shifted towards z < 0 direction (bottom middle panel) and the neck contracts and the warp develops. The warp becomes prominent at t ∼ 5.1 × 10 4 yr. Note that the outflow velocity becomes asymmetric and larger in z > 0 region as the warp develops.
Similar structures have been obtained in some previous studies. Tomida et al. (2013) showed that the warp of pseudo-disk forms in their ideal MHD simulations (although their morphology is asymmetric with respect to z axis). Wurster et al. (2016) and Zhao et al. (2018b) also reported the warps which are very similar to those obtained in our simulations. In turbulent cloud cores, the pseudo-disk tends to have more complicated warp structure as reported by Lam et al. (2019). Note also that Lai (2003) analytically showed that disk-like structure with hour-glass magnetic field can be unstable. Although they investigate the instability of accretion disk and their results can not directly be applied to our configuration, the key mechanism is that the perturbation towards the vertical direction can grow by the Lorentz force. Thus, we speculate the warp structure may be related to this instability. In Appendix B, we show the results of the numerical tests to reinforce the notion that the warp has physical origin.
The warp of the pseudo-disk has a negative impact on disk growth. Due to the warp, the "neck" of the hour-glass magnetic field contracts (see white arrows in figure 13) and magnetic flux in the disk increases. As a result, the magnetic field is strengthened in the disk. Figure 14 shows that magnetic field strength at r ∼ 10 AU increases from 10 −2 G to 10 −1 G during warp development. The magnetic field is vertically density weighted averaged in |z| < 50 AU and is also azimuthally averaged. The magnetic field with warp becomes stronger than that without warp at the corresponding epoch. The increase of the magnetic field strength by the warp enhances the magnetic braking in the disk, and the disk angular momentum and disk size begin to decrease (see figures 15 and 5) after warp formation. Our results show that there is an evolution path in which the disk begins to shrink ∼ 10 4 yr after its birth.
In our simulations, the warp is formed in model MRN, model 01µm, and model a0035µm, and the pseudo-disks in other simulations are approximately symmetric with respect to the x-y plane at least within 10 4 yr after protostar formation. Whether the warp of the pseudo-disk develops in other models in subsequent (long-term) evolution and the detailed physics triggering its formation are unclear and further study is required.

Time evolution of angular momentum and disk radius
In this subsection, we investigate the time evolution of angular momentum and centrifugal radius of the central region to quantitatively discuss disk size evolution after protostar formation.
In this paper, we do not use the disk criteria such as those considering rotation velocity and infall velocity Masson et al. 2016) or those considering gravitational force and centrifugal force (Tsukamoto et al. 2015a) because we found that they overestimate disk size due to contamination from marginally outflowing region when the structure of the inner envelope is elongated due to pseudo-disk warp.
The bottom right panels of figure 2, 4, and 5 highlight this concern. The bottom right panel of figure 5 shows that there is a rapidly rotating low-density region of r ∼ 100 AU on the x-y plane around the central disk. However, from the x-z map of density and plasma β (bottom right panels of figure 2 and figure 4), such low-density and low β regions at r ∼ 100 AU on the x axis do not match with our intuition for the circumstellar disk. Furthermore, as we will show later, the angular momentum of the region decreases after the formation of the warped pseudo-disk, although the disk size estimated by the above-mentioned criteria increases. This strongly suggests that the criteria for disk are inappropriate when the system loses its symmetric structure.
Instead, we use total angular momentum and centrifugal radius of the central region to estimate disk size, which is more physically rigorous and free from disk size overestimation. The angular momentum of disk J(ρ disk ) is calculated as For the density threshold of the disk, we choose ρ disk = 10 −14 g cm −3 . We confirmed that our results below do not strongly depend on the choice of ρ disk (see also, contours of figure 5 and 9). The centrifugal radius is then calculated as  Herej(ρ disk ) = J(ρ disk )/M (ρ disk ) where M (ρ disk ) is the enclosed mass within the region ρ > ρ disk . We regard this centrifugal radius as a disk radius. The solid lines of figure 15 show time evolution of the disk radius. The disk radius at the protostar formation epoch is r ∼ 10 AU and increases in all simulations up to t 5×10 4 yr. However, it begins to decreases in model a0035µm (red), model a01µm (black) and model MRN (yellow). The epochs of disk size decreases corresponds to the epochs of the warp formation. As shown in figure 14, the warp of the pseudodisk strengthens the magnetic field (and magnetic torque) in the disk and has a negative impact on disk growth. In model a03µm and model trMRN, the disk radius as well as J(ρ disk ) continues to increase until the end of the simulation. Among the models, the maximum disk radius of r ∼ 60 AU is realized in model a03µm. The disk size (and angular momentum) evolution of model trMRN are almost identical to those of model a03µm. The evolution of disk radius and angular momentum shows that whether the warp occurs or not affects disk size evolution. Figure 16 shows the time evolution of disk mass M disk (solid) and protostellar mass Mstar (dashed). For disk mass, we used enclosed mass within ρ > ρ disk = 10 −14 g cm −3 . Again, we confirmed that the disk mass does not strongly depends on ρ disk , and ρ disk = 10 −13 g cm −3 and ρ disk = 10 −15 g cm −3 give almost identical results. For the stellar mass, we plot the mass of the sink particle.

Time evolution of mass of protostar and disk
At the protostar formation epoch, the disk mass is M disk ∼ 4 × 10 −2 M , which corresponds to the mass of the pressure-supported first core (Larson 1969;Masunaga et al. 1998), and M disk does not significantly decrease around t = 4.2 × 10 4 yr. Thus, the most part of the first core does not accrete onto the central protostar but stays around the protostar. This means that most of the gas in the first core is directly transformed into the disk. This formation picture of the circumstellar disk from the first core was suggested by  and Inutsuka (2012), and our results are consistent with theirs.
After protostar formation, disk and protostellar mass increase almost monotonically and reach ∼ 0.1 M within 10 4 yr after protostar formation. The slight decrease of disk mass in model a01µm and model MRN is due to enhancement of magnetic braking by the warp of the pseudo-disk. The mass of the disk is comparable to or larger than the protostellar mass within 10 4 yr after protostar formation (or until the central star mass becomes ∼ 0.1 M ). We discuss that a massive disk is a natural consequence of large mass accretion from the envelope in §4.1.
The total mass in the central region is slightly different among the models (∼ 5 × 10 −2 M at t ∼ 5.3 × 10 4 yr, for example). This difference in total mass is consistent with the difference in outflow mass. Mass removal by outflow causes the difference of the total mass in the center.

Time evolution of mass accretion rate
In this subsection, we investigate mass accretion rate, which is a fundamental parameter for the evolution of YSOs. The envelope-to-disk mass accretion rate is calculated aṡ M env,disk = ∆M enclose /∆t interval , where ∆M enclose is the difference of the mass within the region of ρ > 10 −14 g cm −3 (including the mass of the central protostar) during the interval of [t, t + ∆t interval ]. ∆t interval is chosen to be 3 × 10 2 yr. the results were found to be almost unchanged with the other value of ∆t interval and our results discussed below barely depends on the choice of ∆t interval . The diskto-star mass accretion rate is calculated asṀ disk,star = ∆Mstar/∆t interval , where ∆Mstar is the difference in the mass of the protostar during the interval [t, t + ∆t interval ]. Figure 17 shows the mass accretion rates. The solid lines show the mass accretion rate from disk to proto-starṀ disk,star . The mean value during the evolution iṡ M disk,star ∼ 10 −5 M year −1 . This mass accretion rate corresponds to almost the upper limit of the observed mass accretion rate of Class 0 YSOs ). The temporal oscillation ofṀ disk,star is due to mass accretion by the spiral arms. The amplitude of the oscillation is a factor of two to three. The dashed lines show envelope-to-disk mass accretion rate and show thatṀ env,disk ∼ 4 × 10 −5 M year −1 at the protostellar formation epoch and decreases with time. At the end of the simulation, it decreases toṀ env,disk ∼ 10 −5 M year −1 . The oscillation ofṀ env,disk in the latter phase is due to density fluctuation of the outer disk by the spiral arms or warp and whether the oscillation appears or not depends on the choice of the density threshold. Thus, we think the oscillation does not reflect the real change of accretion rate towards the central region. The mass accretion rate ofṀ env,disk ∼ 4 × 10 −5 M year −1 at the protostar formation epoch and its decrease in the latter accretion phase are quantitatively consistent with previous analytic and simulation studies of dynamical collapse of cloud cores (Whitworth & Summers 1985;Saigo & Hanawa 1998;Tomisaka 1996;Vorobyov & Basu 2005). Note that the accretion rate of a dynamical collapse becomes much larger than that of a collapse of a singular isothermal spherė M env,disk ∼ 2 × 10 −6 M year −1 (Shu 1977) because of larger density and larger infall velocity of the envelope.

Time evolution of outflows
In this subsection, we investigate the properties of outflows. We define an outflow as a region that satisfies vr = (v · r)/|r| > 2cs,iso and ρ < 10 −14 g cm −3 where cs,iso = 0.19 km s −1 is sound velocity at T = 10 K. Outflows are formed in all models.
The top left panel of figure 18 shows time evolution of outflow size. The outflow size is defined as the distance of the particle in the outflow that is farthest from the central protostar. In our simulations, strong outflows with velocity of v > 1 km s −1 form several thousand years after protostar formation (note that the protostars form at t = 4.15 × 10 4 yr). Outflow size monotonically increases and reaches ∼ 3000 AU in ∼ 5 × 10 3 yr after outflow formation. This suggests that the mean velocity of the outflow head is ∼ 3 km s −1 .
The top right panel of figure 18 shows the time evolution of the outflow mass Mout. Mout monotonically increases in all models. The outflow mass is anti-correlated with the disk size, suggesting that the outflow activity is related to disk growth The outflow masses and dynamical timescales obtained in our simulations are consistent with observed outflows with dynamical timescale t dyn < 10 4 yr, whose mass ranges from 10 −2 M to 10 −1 M (Wu et al. 2004).
The bottom left panel shows the linear momentum of outflow Pout. The difference in linear momentum mainly comes from the difference of mass of outflow, and the mean velocity of outflow ≡ Pout/Mout is ∼ 2 km s −1 in all simulations. This estimate is consistent with the outflow velocity estimated from the size and the age.
The bottom right panel shows the outflow angular momentum. The outflow angular momentum also shows anticorrelation to the disk size (compare this panel with figure  15) which is consistent with Wurster et al. (2016). The difference in outflow angular momentum among the models is an order of magnitude (e.g., at t ∼ 5.3×10 4 yr) and does not merely come from the difference in mass. The outflow angular momentum is comparable to or even larger than disk angular momentum in t 5 × 10 4 yr for model a0035µm, model a01µm, and model MRN. Coincidently, the disk size begins to decrease at this epoch. Thus, the angular momentum removal due to outflow may play a role in disk size evolution for these models. Note also that pseudo-disk warping occurs in these models. On the other hand, in model a03µm and model trMRN, the outflow angular momentum remains smaller than the disk angular momentum. This indicates that angular momentum removal by outflow plays a minor role for disk evolution in these models.

Magnetic braking
In all models considered in this paper, circumstellar disks are formed immediately after protostar formation. The mass and size of the circumstellar disk are ∼ 10 AU and ∼ 4 × 10 −2 M at its formation epoch, respectively. The initial size and mass are largely consistent with those of the first core (Larson 1969;Masunaga et al. 1998), indicating that the first core directly transforms into the circumstellar disk Inutsuka 2012). The disk grows to several tens of AU at 10 4 yr after protostar formation (protostellar mass reaches M ∼ 0.1 M at this epoch).
The magnetic braking catastrophe (Mellon & Li 2008), which claims that disk formation is completely suppressed by magnetic braking in an early phase of protostar evolution, has been a long-standing issue in theoretical studies on protostar formation (e.g., Allen et al. 2003;Price & Bate 2007b;Mellon & Li 2008;Hennebelle & Fromang 2008;Hennebelle & Ciardi 2009;Joos et al. 2012;Santos-Lima et al. 2012;Seifried et al. 2013;Joos et al. 2013;Li et al. 2013;Krasnopolsky et al. 2011;Li et al. 2011;Tomida et al. 2013Tomida et al. , 2015Tsukamoto et al. 2015b,a;Masson et al. 2016;Wurster et al. 2016;Zhao et al. 2018b;Lam et al. 2019). The key physical mechanisms of magnetic braking catastrophe are magnetic flux freezing and a rapid increase in the magnetic field towards the central star. Magnetic diffusion is hence the most promising mechanism to overcome catastrophic magnetic braking because it relaxes the heart of the magnetic braking catastrophe, i.e., the flux freezing between the magnetic field and gas. Therefore, in principle, the magnetic braking catastrophe can be solved by magnetic diffusion. Our previous studies with three dimensional non-ideal MHD simulations have shown that disk formation (with size of ∼ 1 to 10 AU) actually becomes possible immediately after protostar formation by ohmic and ambipolar diffusion (see the comparison between ideal and resistive MHD simulation in Tsukamoto et al. 2015b). The angular momentum shown in figure 6 of Tsukamoto et al. (2015b) is almost the same as that at the protostar formation epoch shown in figure 15. The simulations of our previous studies were, however, halted just after protostar formation. Therefore, whether the disk can survive and grow after protostar formation was unclear.
The current study shows that the disk grows to several 10 AU scale and is long-lived, at least for ∼ 10 4 yr after protostar formation. Recent theoretical studies considering magnetic diffusion have also confirmed that the disk is formed in a very early phase of protostar formation Wurster et al. 2016;Masson et al. 2016). Furthermore, many observational studies have shown that circumstellar disks form in a very early phase of protostar formation (e.g., Murillo et al. 2013;Ohashi et al. 2014;Yen et al. 2017). Thus, the statement that magnetic braking is catastrophic, in the sense that magnetic braking completely suppresses early disk formation, is not supported either theoretically or observationally.
However, note that the strength of magnetic braking depends on many factors, such as initial conditions, included physics, and microscopic chemistry. Furthermore, the treatment of the central protostar (or inner boundary condition) differs among theoretical studies. This may cause quantitative differences in disk size evolution. Thus, in some cases, a disk did not form even with the non-ideal MHD effect, as indicated in Li et al. (2011) and Zhao et al. (2018b) (see also Machida et al. 2014;Hennebelle et al. 2020, for the impact of the inner boundary condition or sink). This is not surprising because of the differences in the other factors. Note also that magnetic braking certainly has a negative impact on disk growth. Compared to our previous studies on disk formation of unmagnetized cloud core (e.g., Tsukamoto & Machida 2011Tsukamoto et al. 2013bTsukamoto et al. , 2015c, the disk size is small and disk fragmentation is suppressed by the magnetic field. For example, as shown in figure 1 of Tsukamoto & Machida (2011), disk fragmentation is expected with our cloud core (α therm = 0.4 and βrot = 0.03) if the magnetic field is ignored.

Impact of dust size on disk formation and evolution
Recently, it has been suggested that the removal of small dust grains (or dust growth) enhances disk formation (Zhao et al. 2016(Zhao et al. , 2018b. We have confirmed this conclusion. As shown in figure 15, an increase in dust size certainly enhances disk growth. This is due to the enhancement of ambipolar diffusion at ρ ∼ 10 −14 g cm −3 (see figure 1). However, we also found that the difference in dust size does not qualitatively change the disk formation, i.e., it does not determine whether the disk forms or not. Our numerical simulations showed that even with small dust grains, such as in model MRN or model a0035µm, the disk does form immediately after protostar formation and survives.
The occurrence of magnetic field drift in the envelope, on the other hand, is different among the simulations with large and small dust grains. This will change the magnetic field strength of the circumstellar disk by changing the amount of brought-in magnetic flux to the disk. The difference in magnetic field strength in the disk may affect the subsequent long-term evolution of the disk (and possibly MRI activity in the disk). Thus, longer-term simulations ( 10 5 yr after protostar formation) with various dust models would be an important subject for future study.

Disk size evolution and comparison with the observations
Recent observations have revealed that the circumstellar disk is formed in the early evolution phase of YSOs (e.g., Murillo et al. 2013;Ohashi et al. 2014;Aso et al. 2015Aso et al. , 2017Yen et al. 2017). Our results are qualitatively consistent with these results. Then, are the simulation results quantitatively consistent with observations?
To answer this question, we plot the disk size of Class 0/I YSOs from Yen et al. (2017) and disk size evolution obtained in this study in figure 19. The horizontal axis shows the sum of the disk and central star mass because the mass of the protostar of the observations is estimated from the Keplerian rotation velocity at the disk edge or infall velocity, and the contribution of the mass in the disk should also be included. Figure 19 shows that, in the late phase (M > 0.1 M ), our results with large dust grain size are approximately consistent with the observations. For example, model a03µm (orange) and model trMRN (blue) have the almost same disk size of L1527 IRS. On the other hand, the disk size in the simulations with pseudo disk warp is generally smaller than the observational results. However, note that the rotationally supported (marginally outflowing) region extends to 100 AU in these simulations (figure 5) and this region can be observationally regarded as rotationally supported disk, and this possibly explains the discrepancy. In earlier phase (M < 0.1 M ), disk size of the simulation tends to be larger than that of B335. This may be due to the difference of initial angular momentum profile between our initial conditions and the conditions of the real cloud cores. The early evolution phase of the disk is more sensitive to the initial angular momentum profile, and a more realistic velocity field such as turbulence would be suitable to investigate early phase disk evolution (Santos-Lima et al. 2012;Joos et al. 2013;Matsumoto et al. 2017;Lewis & Bate 2018;Lam et al. 2019;Takaishi et al. 2020).

Massive disk formation as a consequence of large mass accretion rate
Our results show that the disks formed in our simulations tends to be massive enough to develop gravitational instability (or Q ∼ 1) where Q is Toomre's Q parameter. Here after "massive disk" is used to mean the marginally gravitationally unstable disk (or disk with Q ∼ 1). It is well known that a massive disk is formed in unmagnetized cloud cores (Nakamoto & Nakagawa 1994;Matsumoto & Hanawa 2003;Vorobyov & Basu 2006;Vorobyov 2009;Vorobyov & Basu 2010;Machida et al. 2010;Tsukamoto & Machida 2011;Stamatellos et al. 2012;Tsukamoto & Machida 2013;Tsukamoto et al. 2013b;Takahashi et al. 2013;Lomax et al. 2014;Tsukamoto et al. 2015c). Our results as well as recent theoretical studies have shown that even with a magnetic field, disk becomes massive and gravitationally unstable once it grows to several 10 AU Tsukamoto et al. 2015a,b;Masson et al. 2016;Zhao et al. 2018b). The large mass accretion rate of 10 −6 M yr −1 causes the formation of the massive disk. This can be understood using the steady viscous accretion disk model as follows. In the viscous accretion disk, the mass accretion rate in the disk, temperature, and surface density, and viscous parameter α are related to where Σgas, cs, Ω, andṀ disk are the gas surface density, sound velocity, orbital period, and mass accretion rate of the circumstellar disk, respectively. This can be rewritten as where Q is Toomre's Q parameter and we approximate epicycle frequency κ is κ = Ω. α is determined by disk temperature and the Q value for a given mass accretion rate of disk. This indicates that even with a massive disk with Q ∼ 1, α ∼ 0.1 is required to realizeṀ disk ∼ 3 × 10 −6 M yr −1 , which is expected value of mass accretion rate in the early evolutionally phase of YSOs (see, Yen et al. 2017). If the disk is less massive (has large Q value), a large value of α is required to realize a mass accretion rate ofṀ disk ∼ 3 × 10 −6 M yr −1 . For example, if Q ∼ 10 in the early disk evolution stage, α as large as 1.2 is required. However, α 1 (meaning that trans-or super-sonic accretion) as well as T 30 K may be unrealistic for a disk with a size of several 10 AU to 100 AU around a low mass protostar. One may think that the mass accretion rate in the disk is not necessarily radially constant and high mass accretion rate only at the inner hot region of disk explains the mass accretion rate of the protostar. However, in Class 0/I YSOs, the gas is continuously supplied from the envelope to the disk with a mass accretion rate ofṀ env,disk ∼ 10 −6 M yr −1 . Thus, ifṀ disk varies in the disk andṀ disk in the outer region of disk is small, the gas stagnates in the outer region and disk mass increases due to envelope accretion. This simple estimate suggests that a massive disk forms when mass accretion rate of 10 −6 to 10 −5 M yr −1 and disk size of several 10 AU to 100 AU are simultaneously realized.
The disk mass estimated from the observations of Class 0 YSOs does not strongly contradict the mass of a marginally gravitationally unstable disk. The disk surface density and disk mass with Q ∼ 2 are estimated as ΣGI(r) = 4.1 Mstar MGI = 2πrΣGI(r)dr ∼ 1.0 × 10 −1 Mstar where we assume T = 150 (r/1AU) −3/7 [K] (Kusaka et al. 1970;Chiang & Goldreich 1997), Keplerian rotation, and Q = const in rin < r < rout. rin and rout are the inner and outer radii of the gravitationally unstable region, re-spectively. We set the Q value for the marginally unstable disk to Qcrit = 2 because the spiral arms develop at Q ∼ 1.4 (Laughlin & Bodenheimer 1994) and the marginally unstable disk may have a slightly larger value than 1.4. Thus, a gravitationally unstable disk with radius of ∼ 100 AU has mass of ∼ 0.1 M . On the other hand, Jørgensen et al. (2009) estimated that the disk mass of Class 0 YSOs has a mean value of ∼ 0.05 M ranging from 0.01 − 0.46 M . Enoch et al. (2011) estimated that the disk mass has a mean value of ∼ 0.2 M ranging from 0.04 − 0.28 M , except for one significantly massive disk. These two studies did not resolve the disks. More recently, Segura-Cox et al. (2018) reported the disk mass of Class 0 YSOs ranging from 0.03 − 0.46 M except for one significantly massive disk in NGC1333IRAS4A (we omit asymmetric objects in the paper). ALMA observations, on the other hand, reported smaller value of disk mass for Class 0/I YSOs. The disk mass of L1527 (Ohashi et al. 2014) and Lupus3 MMS ) are estimated as 1.3 × 10 −2 M and 1.0 × 10 −1 M , respectively. These estimated values (especially L1527 IRS) is smaller than the mass of gravitationally unstable disk. Note, however, that there are several uncertainties in the estimate of disk mass. Dust opacity depends on dust size, composition, and the shape (e.g., Miyake & Nakagawa 1993;Ossenkopf & Henning 1994;Birnstiel et al. 2018). Dust growth and subsequent dust radial drift possibly decreases the dust-to-gas mass ratio and cause underestimation of gas mass (Tsukamoto et al. 2017b). Dust scattering due to the dust growth may also decrease the dust thermal emission, especially for the short wave length (Miyake & Nakagawa 1993;Zhu et al. 2019). Considering these uncertainties, we think that the formation of a massive disk in the early evolution phase does not strongly contradict current observations. However, longer term simulations and detailed comparisons with observations are important subjects for future study.

Condition of drift of magnetic field and ion in the envelope
In §3.3, we investigated magnetic field drift induced by ambipolar diffusion in the envelope. We showed that magnetic field drift more easily occurs with relatively large dust grains in which the ambipolar diffusion is strong in the envelope (see figure 1). We also confirmed that the Hall parameter is β Hall > 1 in r 100 AU in the simulations with magnetic field drift. The Hall parameter indicates the degree of ion-magnetic field coupling and β Hall 1 means that the ion is well coupled to the magnetic field (for details, see Nakano et al. 2002). Thus, our results suggest that ionneutral drift may occur in the envelope of Class 0/I YSOs, especially with relatively large dust grains. Figure 12 shows that the region with outward magnetic field drift has the size of r drift 100 AU when Mstar + M disk 0.1 M and expands to r drift > 100 AU as the mass in the center increases.
Recently, Yen et al. (2018) attempted to observe ionneutral drift in young Class 0 YSO B335. They did not detect ion-neutral drift (more precisely, the velocity difference of ion and neutral is less than δv < 0.3 km s −1 ). Possible explanations for this non-detection is that B335 is too young (Mstar ∼ 0.05 M ) or that the dust size in the enve-lope is not large enough. We suggest that a more evolved Class 0/I YSOs with a central protostar (plus disk) mass of M > 0.1 M may be a good candidate to observe ion-neutral drift. If a velocity difference of ion and neutral is observed in future observations, it will provide a unique opportunity to quantify the magnetic field strength because ion-neutral drift velocity is a function of magnetic field strength.

Formation and early evolution of outflow
In our simulations, outflows are ubiquitously formed. At the end of the simulation, its size reaches several 1000 AU ( figure  18). We found that the mass, linear momentum, and angular momentum of outflow depend on the dust model. The outflows in small dust models (model a0035µm, model a01µm, and model MRN) tend to have larger mass, linear momentum, and angular momentum than those in large dust models (model a03µm and model trMRN). The difference of ηA in the low-density outflow region may cause this difference.
The mass of outflow formed in our simulation is in the range of 10 −2 M < M < 10 −1 M ( figure 18). This is in good agreement with observations. Wu et al. (2004) reported that an observed outflow with dynamical time of 10 4 yr has a mass of 10 −3 M < M < 10 −1 M and mostly within the range of 10 −2 M < M < 10 −1 M . The dispersion of the outflow mass reported in Wu et al. (2004) is possibly due to the difference in ionization degree and hence difference of typical dust size.

SUMMARY
Our results are summarized as follows.
(i) Circumstellar disks are formed in all simulations. The disks have sizes of ∼ 10 AU at the protostar formation epoch and grow to several tens of AU at ∼ 10 4 yr after protostar formation. Disk sizes are almost identical among the simulations as long as pseudo-disk warp does not develop. Once pseudo-disk warp develops, the disk begins to shrink.
(ii) Magnetic field drift in the envelope may occur in the early evolution of young stellar objects. The Hall parameter in the envelope is generally β Hall 1, and ionneutral drift is also expected there. Ion-neutral field drift of v drift 0.19 km s −1 at r 100 AU occurs under conditions with relatively large dust grains of a 0.2µm (or absence of small grain) and protostar (plus disk) mass of M 0.1 M .
(iii) The mass of the circumstellar disk tends to be comparable to the mass of the central star, and gravitational instability develops in the early phase of disk evolution. A massive disk is a consequence of the high mass accretion rate at the early evolution stage.
(iv) The warp of the pseudodisk can develop at ∼ 10 4 yr after protostar formation. The warp enhances magnetic field strength and magnetic braking in the disk, and has a negative impact on disk growth.
(v) Outflows are ubiquitously formed. In some simulations, its angular momentum becomes comparable to the disk angular momentum, and outflow may have a major impact on disk growth.

ACKNOWLEDGMENTS
We thank Dr. Iwasaki Kazunari and Dr. Okuzumi Satoshi for fruitful discussions. We also thank anonymous referee for helpful comments. The computations were performed on a parallel computer, XC40/XC50 system at CfCA of the NAOJ. This work is supported by JSPS KAKENHI grant number 17H06360, 18H05437, 18K13581, 18K03703.

APPENDIX A: INITIAL DENSITY AND MAGNETIC FIELD CONFIGURATION
In this appendix, we describe our initial and boundary conditions in detail, as well as the motivations for adopting the initial conditions.
In long-term simulations of cloud core collapse after protostar formation, the outflows grow to the scale of 10 3 AU i.e., comparable to the initial radius of the core. Thus, precise care is required for the outer boundary. In our previous studies, we adopted a rigidly rotating shell at r ∼ Rc, where Rc is the radius of the cloud core. However, our numerical experiences has shown that such a boundary reflects the outflow and shakes up the density structures in r < Rc, which is clearly numerical. Thus, we need more appropriate outer boundary condition.
Our strategy follows that of , i.e., setting the outer boundary far from the cloud core surface by adding surrounding medium to the core (see e.g., Price & Bate 2007a, for different strategy to impose outer boundary with SPH).  placed the molecular cloud core in a medium with a constant density. The size of medium was 2 5 Rc. With a nested grid code (or AMR code), the outer medium only requires acceptable computational costs. However, with the SPH scheme, the computational cost is proportional to the mass, and hence volume for a constant density medium, and a 2 5 times larger mass requires unacceptably large computational costs.
To avoid this problem, we adopted the Bonner-Evert sphere surrounded by a medium with a steep density profile of ρ ∝ r −4 for r > Rc as described in §2.3. With this profile the total mass of the entire domain is ∼ 2Mc, even when we set the boundary radius to R b = 10Rc.
A problem then arises for the magnetic field structure. If we adopt a constant magnetic field with our density profile, the plasma β obeys β ∝ r −4 in the outer medium, and a low β region emerges, which requires very small time-stepping. To avoid small time-stepping, we constructed a magnetic field profile which has a constant vertical component in the central region, and decreases in the outer region. With this magnetic field profile, plasma β becomes constant in the larger radius and the low β problem is avoided.
The magnetic field of our initial conditions is generated from the vector potential in cylindrical coordinate (R, z) of and the resultant magnetic field is where B0 is the magnetic field at the center. This magnetic field profile has the desired nature. In R → 0, the magnetic field becomes constant and has only the z component. In the spherical coordinate (r, θ), the magnetic field strength is given as ) 2 cos 2θ and except at the midplane (θ = π/2), B ∝ r −2 as r → ∞.
On the other hand, B ∝ r −4 as r → ∞ at the midplane. The ratio of the R and z components of the magnetic field is given as and BR Bz = tan θ = R z (r → ∞).
Thus, the magnetic field is parallel to the position vector (B r) as r → ∞ apart from at the midplane. On the other hand, BR = 0 at the midplane and the magnetic field has only a vertical component at the midplane. With our magnetic field configuration, the magnetic flux is 1/2 times smaller at the edge of the core R = R0 than that with constant magnetic field with B0 although the central magnetic field strength is the same.

APPENDIX B: NUMERICAL TESTS ON THE ORIGIN OF THE WARP
As we have shown, the warp of the pseudodisk often develops in our simulations. To confirm that the warp is not due to the numerical artifact but is physical, we conducted two numerical tests. In this appendix, we describe the results of the numerical tests and show that the warp develops even without a sink particle and even in a simulation with a nested-grid code.
When we first obtained the warp, we were concerned that the numerical artifact of the sink particle possibly causes the warp. To deny this possibility, we conducted the simulation without the sink particle but employing stiff EOS (in other words, we keep using the equation (4) which is not appropriate in high density of ρ 10 −11 g cm −3 ).
The density cross-section of the simulation is shown in figure B1. The initial condition and the dust model are the same as the model MRN. The figure clearly shows that the warp also develops even without a sink particle, although we find that the epoch of warp formation is slightly (∼ 10 3 yr) delayed. Thus, we conclude that the warp is not caused by the numerical artifact of the sink particle.
Another concern was the possible artifact due to the numerical scheme. Although the Godunov SPMHD scheme passes major numerical tests very well (see, Iwasaki & Inutsuka 2011) and we are sure that it can reasonably capture the evolution of the cloud core, the additional test with another numerical scheme was desired.
For this purpose, we conducted a disk formation simula- Figure B1. Density cross-sections on the x-z plane for central 500-AU square region of the simulations with the stiff EOS before and after the development of the warp. The initial conditions and dust model are the same as model MRN.
tion with the numerical code which used in Machida & Basu (2019). The simulation settings of this simulation are different from those of the other simulations presented in this paper. We briefly describe the settings of the simulation below.
The simulation was done with the nested grid code which has been developed by Machida and his collaborators. As the initial condition, the cloud core with a Bonnor-Ebert density profile was adopted. The initial cloud core has a radius of 5.9 × 10 3 AU and a mass of 1 M . A uniform magnetic field of B0 = 43 µG and a rigid rotation of Ω0 = 1.5 × 10 −15 s −1 are added to the initial cloud core, which correspond to the normalized mass-to-flux ratio of µ = 3 and the ratio of rotational to gravitational energy of βrot = 0.02. As the cloud collapses, a finer grid is automatically generated to ensure the Truelove condition, in which the Jeans wave length is resolved at least 16 cells. Each rectangular grid has cells of (i, j, k) = (64, 64, 64). The sink cell technique was used with a threshold density of 10 14 cm −3 and a sink accretion radius 0.5 AU. The grid size and cell width of the finest grid are 24 AU and 0.36 AU, respectively. Equations (1)- (7) in Machida et al. (2020) were solved for this simulation. The figure B2 shows the simulation results and the warp is also formed in this simulation although its size is small compared to the other simulations in this paper, which may reflect the fact that the magnetic field is stronger and the rotation is weaker than the other simulations in this paper. Despite we employed completely different numerical scheme and initial condition, we reproduced the warp. We believe that the reproduction of the warp in these numerical tests strengthens the claim that the warp has physical origin.