Early Evolution of Disk, Outflow, and Magnetic Field of Young Stellar Objects: Impact of Dust Model

In our numerical simulations, non-ideal magnetohydrodynamics (MHD) equations were solved:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{D{\boldsymbol{v}}}{{Dt}} & = & -\displaystyle \frac{1}{\rho }\left\{{\rm{\nabla }}\left(P+\displaystyle \frac{1}{2}| {\boldsymbol{B}}{| }^{2}\right)-{\rm{\nabla }}\cdot ({\boldsymbol{BB}})\right\}\\ & & -{\rm{\nabla }}{\rm{\Phi }},\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{D{\boldsymbol{B}}}{{Dt}} & = & \left({\boldsymbol{B}}\cdot {\rm{\nabla }}{\boldsymbol{v}}-{\boldsymbol{B}}({\rm{\nabla }}\cdot {\boldsymbol{v}})\right)-{\rm{\nabla }}\times \left\{{\eta }_{{\rm{O}}}({\rm{\nabla }}\times {\boldsymbol{B}})\right.\\ & & -\left.{\eta }_{{\rm{A}}}(({\rm{\nabla }}\times {\boldsymbol{B}})\times \hat{{\boldsymbol{B}}})\times \hat{{\boldsymbol{B}}}\right\},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}=4\pi G\rho ,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&P=P(\rho )={c}_{{\rm{s}},\ \mathrm{iso}}^{2}\rho \left\{1+{\left(\displaystyle \frac{\rho }{{\rho }_{\mathrm{crit}}}\right)}^{2/3}\right\}.\end{eqnarray} \tag{ 4 }$$

Here ρ is the gas density, P is the gas pressure, ${\boldsymbol{B}}$ is the magnetic field, and $\hat{{\boldsymbol{B}}}$ is defined as $\hat{{\boldsymbol{B}}}\equiv {\boldsymbol{B}}/| {\boldsymbol{B}}|$. ${\eta }_{{\rm{O}}}$ and ${\eta }_{{\rm{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. ${c}_{{\rm{s}},\mathrm{iso}}=190\,{\rm{m}}\,{{\rm{s}}}^{-1}$ is the isothermal sound velocity at 10 K. We used a critical density of ${\rho }_{\mathrm{crit}}=4\,\times {10}^{-14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ above which gas behaves adiabatically. In this study, we ignored the Hall effect.

We used the smoothed particle magnetohydrodynamics (SPMHD) method to solve the equations (Iwasaki & Inutsuka 2011, 2013). Our numerical code was parallelized with a message-passing interface. We treated ohmic and ambipolar diffusion according to the prescription in Tsukamoto et al. (2013a). The two diffusion processes were accelerated by super-time stepping (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 ${\rho }_{\mathrm{sink}}=4\times {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. In our simulations, one sink particle was permitted. The sink particle absorbs smoothed particle hydrodynamics (SPH) particles with $\rho \gt {\rho }_{\mathrm{sink}}$ within $r\lt {r}_{\mathrm{sink}}=1\,\mathrm{au}$, and the mass and linear momentum of SPH particles are added to those of the sink particle. The sink particle interacts with SPH particles via gravity. The system was integrated up to $\sim {10}^{4}$ yr after protostar formation, at which the mass of the sink particle (or protostar) reached $\sim 0.1{M}_{\odot }$.

We used the tabulated resistivity calculated by the methods described in Susa et al. (2015). We considered the ion species of ${{\rm{H}}}_{3}^{+},{{\rm{H}}}_{2}^{+},{{\rm{H}}}_{3}^{+},{\mathrm{HCO}}^{+},{\mathrm{Mg}}^{+}$, $\ {\mathrm{He}}^{+},{{\rm{C}}}^{+},{{\rm{O}}}^{+},{{\rm{O}}}_{2}^{+},{{\rm{H}}}_{3}{{\rm{O}}}^{+},{\mathrm{OH}}^{+},{{\rm{H}}}_{2}{{\rm{O}}}^{+}$, and the neutral species of ${\rm{H}},{{\rm{H}}}_{2},\ \mathrm{He},\ \mathrm{CO},\ {{\rm{O}}}_{2},\ \mathrm{Mg},{\rm{O}},\ {\rm{C}},\ \mathrm{HCO},\ {{\rm{H}}}_{2}{\rm{O}},\ {\rm{and}}\ \mathrm{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+{\gamma }_{T}{(\rho /{\rho }_{c})}^{({\gamma }_{T}-1)})\,{\rm{K}}$, where ${\gamma }_{T}=7/5$. The dust internal density was fixed to be ${\rho }_{d}=2\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. The cosmic-ray ionization was fixed to be ${\xi }_{\mathrm{CR}}={10}^{-17}\,{{\rm{s}}}^{-1}$.

For the dust models, we considered single-sized dust models and models with a size distribution of $n(a)\propto {a}^{-3.5}$. For the single-sized dust models, we considered three dust sizes of $a=0.035$, $0.1$, and $0.3\,\mu {\rm{m}}$. For the models with size distributions, we considered models with minimum and maximum dust sizes of ${a}_{\min }=0.005\,\mu {\rm{m}}$ and ${a}_{\max }=0.25\,\mu {\rm{m}}$ (MRN distribution) and ${a}_{\min }=0.1\,\mu {\rm{m}}$ and ${a}_{\max }=0.25\,\mu {\rm{m}}$ (truncated MRN distribution, which mimics the size distribution of Zhao et al. 2018b). We assumed a fixed dust-to-gas mass ratio of 0.01.

We adopted a density-enhanced Bonnor–Ebert sphere surrounded by medium with a steep density profile of $\rho \propto {r}^{-4}$ as the initial density profile,

$$\begin{eqnarray}&&\rho (r)=\left\{\begin{array}{l}{\rho }_{0}{\xi }_{\mathrm{BE}}(r/a)\,\mathrm{for}\,r\lt {R}_{c}\\ {\rho }_{0}{\xi }_{\mathrm{BE}}({R}_{c}/a){\left(\tfrac{r}{{R}_{c}}\right)}^{-4}\,\mathrm{for}\,{R}_{c}\lt r\lt 10{R}_{c},\end{array}\right.\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&a={c}_{{\rm{s}},\mathrm{iso}}{\left(\displaystyle \frac{f}{4\pi G{\rho }_{0}}\right)}^{1/2},\end{eqnarray} \tag{ 6 }$$

where ${\xi }_{\mathrm{BE}}$ is the nondimensional density profile of the critical Bonnor–Ebert sphere, f is a numerical factor related to the strength of gravity, and ${R}_{c}=6.45a$ is the radius of the cloud core. f = 1 corresponds to the critical Bonnor–Ebert sphere, and the core with $f\gt 1$ is gravitationally unstable.

A Bonnor–Ebert sphere is determined by specifying central density ${\rho }_{0}$, the ratio of the central density to density at R_c ${\rho }_{0}/\rho ({R}_{c})$, and f. In this study, we adopted the values of ${\rho }_{0}=7.3\times {10}^{-18}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, ${\rho }_{0}/\rho ({R}_{c})=14$, and f = 2.1. Then, the radius of the core is ${R}_{c}=4.8\times {10}^{3}\,\mathrm{au}$, and the enclosed mass within R_c is ${M}_{c}=1{M}_{\odot }$. The ${\alpha }_{\mathrm{therm}}$ ($\equiv {E}_{\mathrm{therm}}/{E}_{\mathrm{grav}}$) is equal to 0.4, where ${E}_{\mathrm{therm}}$ and ${E}_{\mathrm{grav}}$ are the thermal and gravitational energy of the central core (without surrounding medium), respectively. The steep envelope was adopted to place the outer boundary far away from the central cloud core. With the steep profile, the total mass of the entire domain remains $\sim 2{M}_{c}$. For the rotation of the cloud core, we adopted an angular velocity profile of ${\rm{\Omega }}(d)=\tfrac{{{\rm{\Omega }}}_{0}}{\exp (10(d/(1.5{R}_{c}))-1)+1}$, where $d=\sqrt{{x}^{2}+{y}^{2}}$ and ${{\rm{\Omega }}}_{0}=2.3\times {10}^{-13}\,{{\rm{s}}}^{-1}$. ${\rm{\Omega }}(d)$ is almost constant for $d\lt 1.5{R}_{c}$ and rapidly decreases for $d\gt 1.5{R}_{c}$. The ratio of the rotational to gravitational energy ${\beta }_{\mathrm{rot}}$ within the core is ${\beta }_{\mathrm{rot}}$ ($\equiv {E}_{\mathrm{rot}}/{E}_{\mathrm{grav}}$) $=0.03$, where ${E}_{\mathrm{rot}}$ 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⁻² as $r\to \infty$ (see Appendix A for details of the magnetic field configuration). The merit of this profile is that it avoids a low-β region in the surrounding medium, which appears when a constant magnetic field strength and radially decreasing density are adopted. We assumed a characteristic length scale of the field configuration of ${R}_{0}={R}_{c}$. The central magnetic field strength and plasma β were ${B}_{0}=62\,\mu {\rm{G}}$ and $\beta \,=1.6\times {10}^{1}$, respectively. With our magnetic field profile, the mass-to-flux ratio of the core μ relative to the critical value was $\mu /{\mu }_{\mathrm{crit}}=({M}_{c}/{{\rm{\Phi }}}_{\mathrm{mag}})=8$, where ${{\rm{\Phi }}}_{\mathrm{mag}}$ is the magnetic flux of the core and ${\mu }_{\mathrm{crit}}={(0.53/3\pi )(5/G)}^{1/2}$. The mass-to-flux ratio is higher than the observed value because the magnetic field becomes weak in the outer region with our magnetic field configuration. However, note that if we had assumed a constant magnetic field profile with a central magnetic field strength, which has been widely used in previous studies, the mass-to-flux ratio would be ${\mu }_{\mathrm{const}}/{\mu }_{\mathrm{crit}}=4$. This value would be more suitable to compare our magnetic field strength with those of previous studies because we investigated the time evolution of the gas in the central region of the core. We resolved 1 ${M}_{\odot }$ with $3\times {10}^{6}$ SPH particles. Thus, each particle had a mass of $m=3.3\times {10}^{-7}{M}_{\odot }$.

The model names and corresponding dust sizes are summarized in Table 1. Table 1 also summarizes the simulation results.

Table 1.  Model Name, Dust Size, and Summary of the Results

Model name	${a}_{\mathrm{dust}}(\mu {\rm{m}})$	${t}_{\mathrm{star}}(\mathrm{yr})$	${t}_{\mathrm{end}}(\mathrm{yr})$	${r}_{\mathrm{disk},\mathrm{end}}$ (au)	${r}_{\mathrm{outflow},\mathrm{end}}$ (au)	${r}_{\mathrm{drift},\mathrm{end}}$ (au)
model_a0035 $\mu {\rm{m}}$	0.035	$4.15\times {10}^{4}$	$5.6\times {10}^{4}$	29	5100	$\lt 100$
model_a01 $\mu {\rm{m}}$	0.1	$4.15\times {10}^{4}$	$5.3\times {10}^{4}$	18	2300	$\lt 100$
model_a03 $\mu {\rm{m}}$	0.3	$4.15\times {10}^{4}$	5.6 × 104	59	3700	780
model_trMRN	0.1 < a < 0.25	4.15 × 104	5.3 × 104	53	3100	770
model_MRN	0.005 < a < 0.25	4.15 × 104	5.3 × 104	23	2400	< 100

Note. The model with a size distribution has a distribution of $n(a)\propto {a}^{-3.5}$. ${t}_{\mathrm{star}}[\mathrm{yr}]$, ${t}_{\mathrm{end}}[\mathrm{yr}]$, ${r}_{\mathrm{disk},\mathrm{end}}$ [au], ${r}_{\mathrm{outflow},\mathrm{end}}$ [au], and ${r}_{\mathrm{drift},\mathrm{end}}$ [au] are the time at protostar formation, time at the end of the simulation, centrifugal radius at ${t}_{\mathrm{end}}$, outflow size at t_end, and maximum radius at which the magnetic field drift velocity is higher than $0.19\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at t_end, respectively.

Download table as:  ASCIITypeset image

First, we investigate how resistivity depends on the dust model. Figure 1 shows ${\eta }_{{\rm{O}}}$ and ${\eta }_{{\rm{A}}}$ with different dust models. The left panel of Figure 1 shows ${\eta }_{{\rm{O}}}$ and ${\eta }_{{\rm{A}}}$ as a function of density. This indicates that an increase in dust size has two contradicting impacts on ${\eta }_{{\rm{A}}}$. With larger dust size, ${\eta }_{{\rm{A}}}$ is large in the low-density region ($\rho \lesssim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$) and small in the high-density region ($\rho \gtrsim {10}^{-13}\,{\rm{g}}\,{\mathrm{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.

Zoom In Zoom Out Reset image size

The small ${\eta }_{{\rm{A}}}$ with the small dust grains in the low-density 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, the recombination rate in the gas phase decreases. This causes an increase in ion abundance and decrease of ${\eta }_{{\rm{A}}}$ in the low-density region with small dust grains. Note that as dust size increases, ${\eta }_{{\rm{A}}}$ becomes similar to the analytic formula by Shu (1983), ${\eta }_{{\rm{A}}}={B}^{2}/(\gamma C{\rho }^{3/2})$ (dotted line) in the low-density region where $\gamma \,=3.5\times {10}^{13}\,{\mathrm{cm}}^{3}\ {({\rm{g}}{\rm{s}})}^{-1}$ and $C=3\times {10}^{-16}\,{({\rm{g}}{\mathrm{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 the dust-free case. The mutually contradicting impact of dust size on ${\eta }_{{\rm{A}}}$ in the high- and low-density regions may introduce diversity to the evolution of disk, outflow, and magnetic flux.

On the other hand, ${\eta }_{{\rm{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 ${\eta }_{{\rm{A}}}$ as a function of the magnetic field at $\rho =4\times {10}^{-15}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. This shows that the dependence of ${\eta }_{{\rm{A}}}$ on the magnetic field strength is not simple, particularly for small dust sizes. ${\eta }_{{\rm{A}}}$ of a=$0.035\,\mu {\rm{m}}$ is almost independent of the magnetic field in $B\lt {10}^{-2}\,{\rm{G}}$ (and hence behaves like "ohmic diffusion"). As the magnetic field increases, on the other hand, ${\eta }_{{\rm{A}}}$ obeys ${\eta }_{{\rm{A}}}\propto {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 the 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 ${\eta }_{{\rm{A}}}\propto {B}^{2}$ is often assumed for ambipolar diffusion, this is not always valid.

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 are discussed in subsequent sections, are highlighted here.

3.2.1. 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\times {10}^{4}\,\mathrm{yr}$. The top left panel of Figure 2 shows that a weak outflow with $v\lt 1\,\mathrm{km}\,{{\rm{s}}}^{-1}$ is formed. This outflow expands to 150 au in size at $t=4.5\times {10}^{4}\,\mathrm{yr}$ (top middle) but almost stalls (top left). Then a stronger and more collimated outflow with $v\gt 1\,\mathrm{km}\,{{\rm{s}}}^{-1}$ is launched (bottom left). At $t=5.1\times {10}^{4}$ yr (which corresponds to $\sim {10}^{4}\,\mathrm{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 a radius of 1 au. The outflow velocity is roughly determined by rotation velocity (i.e., Keplerian velocity; ${v}_{{\rm{K}}}\sim 1\,\mathrm{km}\,{{\rm{s}}}^{-1}{(100/\mathrm{au})}^{-1/2}{({M}_{\mathrm{star}}/0.1{M}_{\odot })}^{1/2}$ ) at the launching point. This suggests that the fast outflow is launched at several tens of au corresponding to the disk radius.

Zoom In Zoom Out Reset image size

The bottom middle panel of Figure 3 shows that the outflow head reaches ∼1000 au at this epoch. The property of outflow is investigated in more detail in Section 3.8. We find that the warp of the pseudodisk develops at $r\sim 100\,\mathrm{au}$, approximately 10⁴ yr after protostar formation (bottom right panel). This pseudodisk warp strengthens the magnetic field in the disk and negatively impacts disk growth. We revisit the pseudodisk warp in more detail in Section 3.4.

Zoom In Zoom Out Reset image size

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 $\beta \gt {10}^{3}$ in the disk to $\beta \lt {10}^{-2}$ in the upper envelope. At the protostar formation epoch (top left panel), the low-β region is localized in $r\lt 100\,\mathrm{au}$, and it expands as time proceeds.

Zoom In Zoom Out Reset image size

The white arrows around the midplane (around the x-axis) indicate that the magnetic field is highly pinched toward 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 toward $z\lt 0$ and contracts as the warp develops. The warp of the pseudodisk 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 $\beta \gtrsim {10}^{3}$ to $\beta \sim {10}^{2}$ once the warped pseudodisk 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 ($\rho \gtrsim {10}^{-14}\,{\rm{g}}\,{\mathrm{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\sim {10}^{4}\,\mathrm{yr}$ thereafter. Thus, our results are consistent with the disk formation scenario in Machida & Matsumoto (2011) 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 (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 pseudodisk develops, the disk begins to shrink (from the bottom middle to the right panels). A more rigorous analysis of the disk size evolution is presented in Section 3.5.

Zoom In Zoom Out Reset image size

3.2.2. Time Evolution of a Model with Large Dust Grains

Next, we investigate the time evolution of a model with large dust grains, model_a03 $\mu {\rm{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 $\mu {\rm{m}}$, respectively. The clear difference between this model and model_MRN at the protostar formation epoch is the absence of outflow from the 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\gt 1\,\mathrm{km}\,{{\rm{s}}}^{-1}$ eventually forms at $\sim 7\times {10}^{3}\,\mathrm{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 nonaxisymmetric spiral arms in the disk (see Figure 9; a nonaxisymmetric m = 1 mode frequently develops in this model). In this model, the outflow is more collimated and its cylindrical radius ($\sqrt{{x}^{2}+{y}^{2}}$) is smaller than that of model_MRN at $t=5.3\times {10}^{4}\,\mathrm{yr}$. As shown in Section 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 pseudodisk warp. Even though we calculated the system evolution $\sim {10}^{4}$ yr after protostar formation, the pseudodisk 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\times {10}^{4}\,\mathrm{yr}$ (bottom middle panel of Figures 2 and 4). Due to the absence of pseudodisk warp, the large disk is maintained (see Figure 9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 shows the plasma β map of model_a03 $\mu {\rm{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 $\mu {\rm{m}}$ is the thickness of the current sheet at the midplane. The bottom left panel shows that the contours of plasma β at $r\gtrsim 50\,\mathrm{au}$ are sparse around the midplane, indicating that the magnetic field slowly changes toward the vertical direction. Furthermore, the white arrows indicate that the magnetic field is weakly pinched toward the center. On the other hand, the bottom left panel of Figure 4 shows 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 toward the center. This difference comes from the strength of ambipolar diffusion and significant magnetic field drift in the pseudodisk of model_a03 $\mu {\rm{m}}$.

Zoom In Zoom Out Reset image size

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 $\sim {10}^{4}\,\mathrm{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.

Zoom In Zoom Out Reset image size

One of the most important phenomena caused by ambipolar diffusion is the 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 the recent attempt of Yen et al. 2018).

Figure 10 shows the azimuthally averaged radial velocities on the x–y plane at $t=5.3\times {10}^{4}\,\mathrm{yr}$ ($\sim 1.1\times {10}^{4}\,\mathrm{yr}$ after protostar formation). The solid, dashed, and dotted lines show the gas radial velocity ${v}_{{\rm{r}}}$, the radial drift velocity of the magnetic field ${v}_{\mathrm{drfit},{\rm{r}}}$, and the radial velocity of the magnetic field ${v}_{\mathrm{drfit},{\rm{r}}}+{v}_{{\rm{r}}}$, respectively. The drift velocity is calculated as

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{\mathrm{drift}}\equiv {\eta }_{{\rm{A}}}\displaystyle \frac{({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}}{| {\boldsymbol{B}}{| }^{2}}.\end{eqnarray} \tag{ 7 }$$

Among our simulations, model_a03 $\mu {\rm{m}}$ and model_trMRN show significant radial drift in a relatively extended region with a size of $r\gt 100\,\mathrm{au}$. In Figure 10, we also plot ${v}_{\mathrm{dirft},{\rm{r}}}$ of model_MRN as an example of the model with a low drift velocity.

Zoom In Zoom Out Reset image size

In model_a03 $\mu {\rm{m}}$, the drift velocity is positive and reaches $\sim 1\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at $r\sim 100\,\mathrm{au}$ (dashed line). As a result, the total radial velocity of the magnetic field (${v}_{{\rm{r}}}+{v}_{\mathrm{drift},{\rm{r}}}$) is $\sim -0.5\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (dotted line) and is much slower than the gas infall velocity (solid line). We define ${r}_{\mathrm{drift}}$ as the maximum value of radius at which the radial magnetic field drift velocity is higher than $0.19\,\mathrm{km}\,{{\rm{s}}}^{-1}$ on the x–y plane. Here we use a velocity threshold of $0.19\,\mathrm{km}\,{{\rm{s}}}^{-1}$, which corresponds to the sound velocity of T = 10 K and approximates the sound velocity in the envelope. In model_a03 $\mu {\rm{m}}$, ${r}_{\mathrm{drift}}$ is ∼700 au at $t=5.3\times {10}^{4}\,\mathrm{yr}$.

The radial drift velocity becomes low as the typical dust size decreases. In model_trMRN, the drift velocity is typically $\sim 0.4\,\mathrm{km}\,{{\rm{s}}}^{-1}$. On the other hand, ${r}_{\mathrm{drift}}$ also extends to ∼700 au in this model. In model_MRN, no notable outward radial drift can be observed, and ${v}_{{\rm{r}}}+{v}_{\mathrm{drift},{\rm{r}}}$ (dotted line) is almost identical to the gas infall velocity (${v}_{{\rm{r}}}$: solid line). This clearly shows that the outward radial drift of the magnetic field occurs only with relatively large dust grains (or in the absence of small dust grains) in the early evolution phase of YSOs. This is mainly because of the difference of ${\eta }_{{\rm{A}}}$ in the low-density region of $\rho \lesssim {10}^{-14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$.

One of the most important aspects of the magnetic field drift by ambipolar diffusion is that the magnetic field drift accompanies the 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 the drift velocity of an ion (${v}_{\mathrm{drift},\mathrm{ion}}$) and the magnetic field (${v}_{\mathrm{drift}}$) can be calculated as

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{\mathrm{drift},\mathrm{ion}}}{{v}_{\mathrm{drift}}}=\displaystyle \frac{{\beta }_{\mathrm{Hall}}^{2}}{(1+{\beta }_{\mathrm{Hall}}^{2})}(1+\xi {\beta }_{\mathrm{Hall}}^{-1}),\end{eqnarray} \tag{ 8 }$$

where ξ is a constant of order unity, and ${\beta }_{\mathrm{Hall}}$ is the Hall parameter,

$$\begin{eqnarray}&&{\beta }_{\mathrm{Hall}}\equiv \displaystyle \frac{{q}_{i}B}{{m}_{i}c}\displaystyle \frac{{m}_{i}+{m}_{n}}{\rho \langle \sigma v{\rangle }_{i}},\end{eqnarray} \tag{ 9 }$$

where q_i and m_i are the charge and mass of the ion, respectively. m_n is the mean neutral mass. $\langle \sigma v{\rangle }_{i}$ is the rate coefficient for collisional momentum transfer between ion and neutral. For the ion species with ${\beta }_{\mathrm{Hall}}\gg 1$, the ratio of velocity becomes ${v}_{\mathrm{drift},\mathrm{ion}}/{v}_{\mathrm{drift}}\sim 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 $\mu {\rm{m}}$ and of model_trMRN at $t=5.3\times {10}^{4}\,\mathrm{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 that $\langle \sigma v{\rangle }_{i}$ is constant and $\langle \sigma v{\rangle }_{i}=1.8\times {10}^{-9}\,{\mathrm{cm}}^{3}\ \,{{\rm{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 ${\beta }_{\mathrm{Hall}}\gt 1$ in the region of $r\gtrsim 100\,\mathrm{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\gtrsim 100\,\mathrm{au}$.

Zoom In Zoom Out Reset image size

The size of the region in which ion-neutral drift occurs is an important quantity to observe ion-neutral drift. Figure 12 shows the time evolution of ${r}_{\mathrm{drift}}$ as a function of ${M}_{\mathrm{star}}+{M}_{\mathrm{disk}}$ for model_a03 $\mu {\rm{m}}$ and for model_trMRN. Here ${M}_{\mathrm{star}}$ is the mass of the sink particle and ${M}_{\mathrm{disk}}$ is the disk mass, which is defined as the enclosed mass of the region with $\rho \,\gt {10}^{-14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ (see Section 3.6). The figure shows that ${r}_{\mathrm{drift}}$ monotonically increases as the mass (and hence magnetic flux) accumulates to the central region and reaches ∼100 au at ${M}_{\mathrm{star}}+{M}_{\mathrm{disk}}\sim 0.1{M}_{\odot }$. This indicates that ${M}_{\mathrm{star}}+{M}_{\mathrm{disk}}\,\gtrsim 0.1{M}_{\odot }$ is required for magnetic field drift to occur in $r\gt 100\,\mathrm{au}$ and indicates that a sufficient amount of mass (and hence magnetic flux) should have accreted onto 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 Section 4.2.

Zoom In Zoom Out Reset image size

As shown in Figures 2 and 4, an interesting structure develops in the pseudodisk, i.e., the warp of the pseudodisk. In our simulations, the warp appears in model_MRN, model_a01 $\mu {\rm{m}}$ and model_a0035 $\mu {\rm{m}}$ at $\sim {10}^{4}\,\mathrm{yr}$ after protostar formation.

Figure 13 shows the time evolution of the density for model_a01 $\mu {\rm{m}}$ from $t=4.85\times {10}^{4}$ yr (just before the warp develops) to $t=5.1\times {10}^{4}$ yr and the development of the warp. At $t\sim 4.85\times {10}^{4}\,\mathrm{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 the $z\lt 0$ direction (top right panel). The gas accretion flow is also shifted toward the $z\lt 0$ direction (bottom middle panel), and the neck contracts and the warp develops. The warp becomes prominent at $t\sim 5.1\times {10}^{4}\,\mathrm{yr}$. Note that the outflow velocity becomes asymmetric and higher in the $z\gt 0$ region as the warp develops.

Zoom In Zoom Out Reset image size

Similar structures have been obtained in some previous studies. Tomida et al. (2013) showed that the warp in the pseudodisk forms in their ideal magnetohydrodynamics (MHD) simulations (although their morphology is asymmetric with respect to the z-axis). Wurster et al. (2016) and Zhao et al. (2018b) also reported warps that are very similar to those obtained in our simulations. In turbulent cloud cores, the pseudodisk tends to have a more complicated warp structure, as reported by Lam et al. (2019). Note also that Lai (2003) analytically showed that a disk-like structure with an hour-glass magnetic field can be unstable. Although they investigated the instability of an accretion disk and their results cannot directly be applied to our configuration, the key mechanism is that the perturbation toward the vertical direction can grow through the Lorentz force. Thus, we speculate that 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 a phyiscal origin.

The warp of the pseudodisk has a negative impact on disk growth. Due to the warp, the neck of the hour-glass magnetic field contracts (see the white arrows in Figure 13) and the 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\sim 10$ au increases from ${10}^{-2}$ to 10⁻¹ G during warp development. The magnetic field is vertically density weighted averaged in $| z| \lt 50\,\mathrm{au}$ and is also azimuthally averaged. The magnetic field with the warp becomes stronger than thae field without warp at the corresponding epoch. The increase in 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 $\sim {10}^{4}\,\mathrm{yr}$ after its birth.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In our simulations, the warp is formed in model_MRN, model_01 $\mu {\rm{m}}$, and model_a0035 $\mu {\rm{m}}$, and the pseudodisks in other simulations are approximately symmetric with respect to the x–y plane, at least within $\lesssim {10}^{4}\,\mathrm{yr}$ after protostar formation. Whether the warp of the pseudodisk develops in other models in subsequent (long-term) evolution is unclear, as is the detailed physics triggering its formation. Further study is required.

In this subsection, we investigate the time evolution of the angular momentum and the centrifugal radius of the central region to quantitatively discuss the disk size evolution after protostar formation.

In this paper, we do not use disk criteria such as those considering rotation velocity and infall velocity (Machida et al. 2011; Masson et al. 2016) or those considering gravitational force and centrifugal force (Tsukamoto et al. 2015a) because we found that they overestimate the disk size due to contamination from the marginally outflowing region when the structure of the inner envelope is elongated due to pseudodisk warp.

The bottom right panels of Figures 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\sim 100\,\mathrm{au}$ on the x–y plane around the central disk. However, from the x–z map of density and plasma β (bottom right panels of Figures 2 and 4), such low-density and low-β regions at $r\sim 100\,\mathrm{au}$ on the x-axis do not match our intuition for the circumstellar disk. Furthermore, as we show below, the angular momentum of the region decreases after the formation of the warped pseudodisk, although the disk size estimated by the above-mentioned criteria increases. This strongly suggests that the criteria for the disk are inappropriate when the system loses its symmetric structure.

Instead, we use the total angular momentum and the centrifugal radius of the central region to estimate the disk size, which is more physically rigorous and free from disk size overestimation. The angular momentum of the disk $J({\rho }_{\mathrm{disk}})$ is calculated as

$$\begin{eqnarray}&&J({\rho }_{\mathrm{disk}})\equiv \left|{\int }_{\rho \gt {\rho }_{\mathrm{disk}}}\rho ({\boldsymbol{r}}\times {\boldsymbol{v}}){dV}\right|.\end{eqnarray} \tag{ 10 }$$

For the density threshold of the disk, we choose ${\rho }_{\mathrm{disk}}={10}^{-14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. We confirmed that our results below do not strongly depend on the choice of ${\rho }_{\mathrm{disk}}$ (see also the contours in Figures 5 and 9). The centrifugal radius is then calculated as

$$\begin{eqnarray}&&{r}_{\mathrm{disk}}\equiv {r}_{\mathrm{cent}}=\displaystyle \frac{\bar{j}{({\rho }_{\mathrm{disk}})}^{2}}{{{GM}}_{\mathrm{star}}}.\end{eqnarray} \tag{ 11 }$$

Here $\bar{j}({\rho }_{\mathrm{disk}})=J({\rho }_{\mathrm{disk}})/M({\rho }_{\mathrm{disk}})$, where $M({\rho }_{\mathrm{disk}})$ is the mass enclosed within the region $\rho \gt {\rho }_{\mathrm{disk}}$. We regard this centrifugal radius as a disk radius.

The solid lines of Figure 15 show the time evolution of the disk radius. The disk radius at the protostar formation epoch is $r\sim 10\,\mathrm{au}$ and increases in all simulations up to $t\lesssim 5\,\times {10}^{4}\,\mathrm{yr}$. However, it begins to decrease in model_a0035 $\mu {\rm{m}}$ (red), model_a01 $\mu {\rm{m}}$ (black), and model_MRN (yellow). The epochs of disk size decreases correspond 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 $\mu {\rm{m}}$ and model_trMRN, the disk radius as well as $J({\rho }_{\mathrm{disk}})$ continue to increase until the end of the simulation. In the models, the maximum disk radius of $r\sim 60\,\mathrm{au}$ is realized in model_a03 $\mu {\rm{m}}$. The disk size (and angular momentum) evolution of model_trMRN is almost identical to that in model_a03 $\mu {\rm{m}}$. The evolution of disk radius and angular momentum shows that the disk size evolution is affected by whether or not the warp occurs.

Figure 16 shows the time evolution of the disk mass ${M}_{\mathrm{disk}}$ (solid) and protostellar mass ${M}_{\mathrm{star}}$ (dashed). For the disk mass, we used the mass enclosed within $\rho \gt {\rho }_{\mathrm{disk}}={10}^{-14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. Again, we confirmed that the disk mass does not strongly depend on ${\rho }_{\mathrm{disk}}$, and that ${\rho }_{\mathrm{disk}}={10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and ${\rho }_{\mathrm{disk}}\,={10}^{-15}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ give almost identical results. For the stellar mass, we plot the mass of the sink particle.

Zoom In Zoom Out Reset image size

At the protostar formation epoch, the disk mass is ${M}_{\mathrm{disk}}\sim 4\times {10}^{-2}{M}_{\odot }$, which corresponds to the mass of the pressure-supported first core (Larson 1969; Masunaga et al. 1998), and ${M}_{\mathrm{disk}}$ does not significantly decrease around $t=4.2\times {10}^{4}\,\mathrm{yr}$. Thus, the greater 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 Machida & Matsumoto (2011) and Inutsuka (2012), and our results are consistent with theirs.

After protostar formation, disk and protostellar mass increase almost monotonically and reach $\sim 0.1{M}_{\odot }$ within 10⁴ yr after protostar formation. The slight decrease in disk mass in model_a01 $\mu {\rm{m}}$ and model_MRN is due to the enhancement in magnetic braking by the warp of the pseudodisk. The mass of the disk is comparable to or higher than the protostellar mass within 10⁴ yr after protostar formation (or until the central star mass becomes $\sim 0.1{M}_{\odot }$). We discuss that a massive disk is a natural consequence of large mass accretion from the envelope in Section 4.1.

The total mass in the central region is slightly different in the models ($\sim 5\times {10}^{-2}{M}_{\odot }$ at $t\sim 5.3\times {10}^{4}\,\mathrm{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.

In this subsection, we investigate the mass-accretion rate, which is a fundamental parameter for the evolution of YSOs. The rate of the envelope-to-disk mass accretion is calculated as ${\dot{M}}_{\mathrm{env},\mathrm{disk}}={\rm{\Delta }}{M}_{\mathrm{enclose}}/{\rm{\Delta }}{t}_{\mathrm{interval}}$, where ${\rm{\Delta }}{M}_{\mathrm{enclose}}$ is the difference in mass within the region of $\rho \gt {10}^{-14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ (including the mass of the central protostar) during the interval of $[t,t+{\rm{\Delta }}{t}_{\mathrm{interval}}]$. ${\rm{\Delta }}{t}_{\mathrm{interval}}$ is chosen to be $3\times {10}^{2}$ yr. The results were found to be almost unchanged with the other value of ${\rm{\Delta }}{t}_{\mathrm{interval}}$, and our results discussed below barely depend on the choice of ${\rm{\Delta }}{t}_{\mathrm{interval}}$. The rate of disk-to-star mass accretion is calculated as ${\dot{M}}_{\mathrm{disk},\mathrm{star}}={\rm{\Delta }}{M}_{\mathrm{star}}/{\rm{\Delta }}{t}_{\mathrm{interval}}$, where ${\rm{\Delta }}{M}_{\mathrm{star}}$ is the difference in the mass of the protostar during the interval $[t,t+{\rm{\Delta }}{t}_{\mathrm{interval}}]$.

Figure 17 shows the mass-accretion rates. The solid lines show the mass-accretion rate from the disk to the protostar ${\dot{M}}_{\mathrm{disk},\mathrm{star}}$. The mean value during the evolution is ${\dot{M}}_{\mathrm{disk},\mathrm{star}}\sim {10}^{-5}{M}_{\odot }{\mathrm{yr}}^{-1}$. This mass-accretion rate corresponds to almost the upper limit of the observed mass-accretion rate of Class 0 YSOs (Yen et al. 2017). The temporal oscillation of ${\dot{M}}_{\mathrm{disk},\mathrm{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 the rate of envelope-to-disk mass accretion, and they show that ${\dot{M}}_{\mathrm{env},\mathrm{disk}}\sim 4\times {10}^{-5}{M}_{\odot }{\mathrm{yr}}^{-1}$ at the protostellar formation epoch and decreases with time. At the end of the simulation, it decreases to ${\dot{M}}_{\mathrm{env},\mathrm{disk}}\sim {10}^{-5}{M}_{\odot }{\mathrm{yr}}^{-1}$. The oscillation of ${\dot{M}}_{\mathrm{env},\mathrm{disk}}$ in the later 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 in accretion rate toward the central region. The mass-accretion rate of ${\dot{M}}_{\mathrm{env},\mathrm{disk}}\sim 4\times {10}^{-5}{M}_{\odot }{\mathrm{yr}}^{-1}$ at the protostar formation epoch and its decrease in the later accretion phase are quantitatively consistent with previous analytic and simulation studies of the dynamical collapse of cloud cores (Whitworth & Summers 1985; Tomisaka 1996; Saigo & Hanawa 1998; Vorobyov & Basu 2005). Note that the accretion rate of a dynamical collapse becomes much higher than that of a collapse of a singular isothermal sphere ${\dot{M}}_{\mathrm{env},\mathrm{disk}}\sim 2\times {10}^{-6}{M}_{\odot }{\mathrm{yr}}^{-1}$ (Shu 1977) because of the higher density and higher infall velocity of the envelope.

Zoom In Zoom Out Reset image size

In this subsection, we investigate the properties of the outflows. We define an outflow as a region that satisfies ${v}_{{\rm{r}}}=({\boldsymbol{v}}\cdot {\boldsymbol{r}})/| {\boldsymbol{r}}| \gt 2{c}_{{\rm{s}},\mathrm{iso}}$ and $\rho \lt {10}^{-14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, where ${c}_{{\rm{s}},\mathrm{iso}}\,=0.19\,\mathrm{km}\,{{\rm{s}}}^{-1}$ is the sound velocity at T = 10 K. Outflows are formed in all models.

The top left panel of Figure 18 shows the time evolution of the 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 a velocity of $v\gt 1\,\mathrm{km}\,{{\rm{s}}}^{-1}$ form several thousand years after protostar formation (note that the protostars form at $t=4.15\times {10}^{4}$ yr). The outflow size monotonically increases and reaches ∼3000 au in $\sim 5\times {10}^{3}\,\mathrm{yr}$ after outflow formation. This suggests that the mean velocity of the outflow head is $\sim 3\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

The top right panel of Figure 18 shows the time evolution of the outflow mass ${M}_{\mathrm{out}}$. ${M}_{\mathrm{out}}$ monotonically increases in all models. The outflow mass is anticorrelated 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 a dynamical timescale ${t}_{\mathrm{dyn}}\lt {10}^{4}\,\mathrm{yr}$, whose mass ranges from ${10}^{-2}{M}_{\odot }$ to ${10}^{-1}{M}_{\odot }$ (Wu et al. 2004).

The bottom left panel shows the linear momentum of the outflow ${P}_{\mathrm{out}}$. The difference in linear momentum mainly comes from the difference in mass of the outflow, and the mean velocity of the outflow $\equiv {P}_{\mathrm{out}}/{M}_{\mathrm{out}}$ is $\sim 2\,\mathrm{km}\,{{\rm{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 an 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 in the models is one order of magnitude (e.g., at $t\sim 5.3\times {10}^{4}$ yr) and does not merely come from the difference in mass. The outflow angular momentum is comparable to or even larger than the disk angular momentum in $t\gtrsim 5\times {10}^{4}\,\mathrm{yr}$ for model_a0035 $\mu {\rm{m}}$, model_a01 $\mu {\rm{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 the disk size evolution for these models. Note also that pseudodisk warping occurs in these models. On the other hand, in model_a03 $\mu {\rm{m}}$ and model_trMRN, the outflow angular momentum remains smaller than the disk angular momentum. This indicates that angular momentum removal by the outflow plays a minor role for the disk evolution in these models.

4.1.1. Magnetic Braking

4.1.2. Impact of the Dust Size on Disk Formation and Evolution

4.1.3. Disk Size Evolution and Comparison with 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. 2015, 2017; Yen et al. 2017). Our results are qualitatively consistent with these results. The question then arises whether the simulation results are quantitatively consistent with observations.

To answer this question, we plot the disk size of Class 0/I YSOs from Yen et al. (2017) and the 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\gt 0.1{M}_{\odot }$), our results with a large dust grain size are approximately consistent with observations. For example, model_a03 $\mu {\rm{m}}$ (orange) and model_trMRN (blue) have almost the same disk size as L1527 IRS. On the other hand, the disk size in the simulations with pseudodisk 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 a rotationally supported disk. This possibly explains the discrepancy. In an earlier phase ($M\lt 0.1{M}_{\odot }$), the disk size of the simulation tends to be larger than that of B335. This may be due to the difference in 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).

Zoom In Zoom Out Reset image size

4.1.4. Massive Disk Formation as a Consequence of a High Mass-accretion Rate

Our results show that the disks formed in our simulations tends to be massive enough to develop gravitational instability (or $Q\sim 1$), where Q is Toomre's Q parameter. Hereafter, "massive disk" is used to mean the marginally gravitationally unstable disk (or the disk with $Q\sim 1$). It is well known that a massive disk is formed in unmagnetized cloud cores (Nakamoto & Nakagawa 1994; Matsumoto & Hanawa 2003; Vorobyov & Basu 2006, 2010; Vorobyov 2009; Machida et al. 2010; Tsukamoto & Machida 2011, 2013; Stamatellos et al. 2012; Tsukamoto et al. 2013b, 2015c; Takahashi et al. 2013; Lomax et al. 2014). Our results as well as recent theoretical studies have shown that even with a magnetic field, the disk becomes massive and gravitationally unstable once it grows to several 10 au (Machida et al. 2011; Tsukamoto et al. 2015a, 2015b; Masson et al. 2016; Zhao et al. 2018b).

The high mass-accretion rate of $\gtrsim {10}^{-6}\,{M}_{\odot }\,{\mathrm{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 the viscous parameter α are related to

$$\begin{eqnarray}&&\alpha {{\rm{\Sigma }}}_{\mathrm{gas}}\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{\rm{\Omega }}}=\displaystyle \frac{{\dot{M}}_{\mathrm{disk}}}{3\pi },\end{eqnarray} \tag{ 12 }$$

where ${{\rm{\Sigma }}}_{\mathrm{gas}}$, ${c}_{{\rm{s}}}$, Ω, and ${\dot{M}}_{\mathrm{disk}}$ are the gas surface density, the sound velocity, the orbital period, and the mass-accretion rate of the circumstellar disk, respectively. This can be rewritten as

$$\begin{eqnarray}\begin{array}{rcl}\alpha & = & \displaystyle \frac{{\dot{M}}_{\mathrm{disk}}}{3{c}_{{\rm{s}}}^{3}}\displaystyle \frac{{c}_{{\rm{s}}}{\rm{\Omega }}}{\pi {{\rm{\Sigma }}}_{\mathrm{disk}}}=\displaystyle \frac{1}{3}\displaystyle \frac{{\dot{M}}_{\mathrm{disk}}}{{c}_{{\rm{s}}}^{3}/G}Q\\ & = & 0.24\left(\displaystyle \frac{{\dot{M}}_{\mathrm{disk}}}{3\times {10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\left(\displaystyle \frac{Q}{2}\right){\left(\displaystyle \frac{T}{30{\rm{K}}}\right)}^{-3/2},\end{array}\end{eqnarray} \tag{ 13 }$$

where Q is Toomre's Q parameter, and we approximate the epicycle frequency κ as $\kappa ={\rm{\Omega }}$. α is determined by the disk temperature, and the Q value for a given mass-accretion rate of disk. This indicates that even with a massive disk with $Q\sim 1$, $\alpha \sim 0.1$ is required to realize ${\dot{M}}_{\mathrm{disk}}\sim 3\times {10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, which is the expected value of the mass-accretion rate in the early evolutionally phase of YSOs (see, Yen et al. 2017). If the disk is less massive (has a high Q value), a high value of α is required to realize a mass-accretion rate of ${\dot{M}}_{\mathrm{disk}}\sim 3\times {10}^{-6}{M}_{\odot }{\mathrm{yr}}^{-1}$. For example, if $Q\sim 10$ in the early disk evolution stage, α as large as 1.2 is required. However, $\alpha \gtrsim 1$ (meaning trans- or supersonic accretion) as well as $T\gg 30$ K may be unrealistic for a disk with a size of several 10–100 au around a low-mass protostar. One may think that the mass-accretion rate in the disk is not necessarily radially constant and that high mass-accretion rate only in the inner hot region of the 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 ${\dot{M}}_{\mathrm{env},\mathrm{disk}}\sim {10}^{-6}{M}_{\odot }{\mathrm{yr}}^{-1}$. Thus, if ${\dot{M}}_{\mathrm{disk}}$ varies in the disk and ${\dot{M}}_{\mathrm{disk}}$ in the outer region of disk is small, the gas stagnates in the outer region and the disk mass increases due to envelope accretion. This simple estimate suggests that a massive disk forms for a mass-accretion rate of 10⁻⁶–${10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and a disk size of several 10–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\sim 2$ are estimated as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Sigma }}}_{\mathrm{GI}}(r) & = & 4.1{\left(\displaystyle \frac{{M}_{\mathrm{star}}}{0.1{M}_{\odot }}\right)}^{1/2}\\ & & \times {\left(\displaystyle \frac{Q}{2}\right)}^{-1}{\left(\displaystyle \frac{r}{100\mathrm{au}}\right)}^{-12/7}{\rm{g}}\,{\mathrm{cm}}^{-2},\end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{GI}} & = & \displaystyle \int 2\pi r{{\rm{\Sigma }}}_{\mathrm{GI}}(r){dr}\sim 1.0\times {10}^{-1}{\left(\displaystyle \frac{{M}_{\mathrm{star}}}{0.1{M}_{\odot }}\right)}^{1/2}\\ & \times & {\left(\displaystyle \frac{Q}{2}\right)}^{-1}{\left(\displaystyle \frac{{r}_{\mathrm{out}}}{100\mathrm{au}}\right)}^{2/7}\left(1-{\left(\displaystyle \frac{{r}_{\mathrm{in}}}{{r}_{\mathrm{out}}}\right)}^{2/7}\right),{M}_{\odot },\end{array}\end{eqnarray} \tag{ 15 }$$

where we assume $T=150\,{(r/1\mathrm{au})}^{-3/7}[{\rm{K}}]$ (Kusaka et al. 1970; Chiang & Goldreich 1997), Keplerian rotation, and $Q=\mathrm{const}$ in ${r}_{\mathrm{in}}\lt r\lt {r}_{\mathrm{out}}$. ${r}_{\mathrm{in}}$ and ${r}_{\mathrm{out}}$ are the inner and outer radii of the gravitationally unstable region, respectively. We set the Q value for the marginally unstable disk to ${Q}_{\mathrm{crit}}=2$ because the spiral arms develop at $Q\sim 1.4$ (Laughlin & Bodenheimer 1994) and the marginally unstable disk may have a slightly higher value than 1.4. Thus, a gravitationally unstable disk with a radius of ∼100 au has a mass of $\sim 0.1{M}_{\odot }$.

On the other hand, Jørgensen et al. (2009) estimated that the disk mass of Class 0 YSOs has a mean value of $\sim 0.05{M}_{\odot }$ ranging from $0.01\,\mathrm{to}\,0.46{M}_{\odot }$. Enoch et al. (2011) estimated that the disk mass has a mean value of $\sim 0.2{M}_{\odot }$ ranging from $0.04\,\mathrm{to}\,0.28{M}_{\odot }$, except for one significantly massive disk. These two studies did not resolve the disks. More recently, Segura-Cox et al. (2018) reported a disk mass of Class 0 YSOs ranging from $0.03\,\mathrm{to}\,0.46{M}_{\odot }$ except for one significantly massive disk in NGC 1333IRAS4A (we omit asymmetric objects in the paper). Atacama Large Millimeter/submillimeter Array observations, on the other hand, reported a lower value of the disk mass for Class 0/I YSOs. The disk masses of L1527 (Ohashi et al. 2014) and Lupus3 MMS (Yen et al. 2017) are estimated as $1.3\times {10}^{-2}{M}_{\odot }$ and $1.0\times {10}^{-1}{M}_{\odot }$, respectively. These estimated values (especially L1527 IRS) are lower than the mass of a gravitationally unstable disk. Note, however, that there are several uncertainties in the estimate of the 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 decrease the dust-to-gas mass ratio and cause an underestimation of the gas mass (Tsukamoto et al. 2017b). Dust scattering due to dust growth may also decrease the dust thermal emission, especially for the short wavelength (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.

In Section 3.3, we investigated the magnetic field drift induced by ambipolar diffusion in the envelope. We showed that the 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 ${\beta }_{\mathrm{Hall}}\gt 1$ in $r\gtrsim 100\,\mathrm{au}$ in simulations with a magnetic field drift. The Hall parameter indicates the degree of ion-magnetic field coupling, and ${\beta }_{\mathrm{Hall}}\gg 1$ means that the ion is well coupled to the magnetic field (for details, see Nakano et al. 2002). Thus, our results suggest that an ion-neutral drift may occur in the envelope of Class 0/I YSOs, especially with relatively large dust grains. Figure 12 shows that the region with an outward magnetic field drift has a size of ${r}_{\mathrm{drift}}\lesssim 100\,\mathrm{au}$ when ${M}_{\mathrm{star}}+{M}_{\mathrm{disk}}\lesssim 0.1{M}_{\odot }$ and expands to ${r}_{\mathrm{drift}}\gt 100$ au as the mass in the center increases.

Recently, Yen et al. (2018) attempted to observe the ion-neutral drift in young Class 0 YSO B335. They did not detect the ion-neutral drift (more precisely, the velocity difference of ions and neutrals is lower than $\delta v\lt 0.3\,\mathrm{km}\,{{\rm{s}}}^{-1}$). Possible explanations for this nondetection are that B335 is too young (${M}_{\mathrm{star}}\sim 0.05{M}_{\odot }$) or that the dust size in the envelope is not large enough. We suggest that a more evolved Class 0/I YSOs with a central protostar (plus disk) mass of $M\gt 0.1{M}_{\odot }$ may be a good candidate to observe the ion-neutral drift. If a velocity difference of ions and neutrals is observed in future observations, it will provide a unique opportunity to quantify the magnetic field strength because the ion-neutral drift velocity is a function of the magnetic field strength.

In our simulations, outflows are ubiquitously formed. At the end of the simulation, their sizes reach 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 $\mu {\rm{m}}$, model_a01 $\mu {\rm{m}}$, and model_MRN) tend to have higher mass, linear momentum, and angular momentum than those in large-dust models (model_a03 $\mu {\rm{m}}$ and model_trMRN). The difference of ${\eta }_{{\rm{A}}}$ in the low-density outflow region may cause this difference.

The mass of the outflow that formed in our simulation is in the range of ${10}^{-2}{M}_{\odot }\lt M\lt {10}^{-1}{M}_{\odot }$ (Figure 18). This is in good agreement with observations. Wu et al. (2004) reported that an observed outflow with a dynamical time of $\lesssim {10}^{4}\,\mathrm{yr}$ has a mass of ${10}^{-3}{M}_{\odot }\lt M\lt {10}^{-1}{M}_{\odot }$ and lies mostly within the range of ${10}^{-2}{M}_{\odot }\lt M\lt {10}^{-1}{M}_{\odot }$. The dispersion of the outflow mass reported in Wu et al. (2004) is possibly due to the difference in ionization degree and hence to the difference in typical dust size.