Radiation Magnetohydrodynamic Simulations of Protostellar Collapse: NonIdeal Magnetohydrodynamic Effects and Early Formation of Circumstellar Disks
Abstract
The transport of angular momentum by magnetic fields is a crucial physical process in formation and evolution of stars and disks. Because the ionization degree in star forming clouds is extremely low, nonideal magnetohydrodynamic (MHD) effects such as ambipolar diffusion and Ohmic dissipation work strongly during protostellar collapse. These effects have significant impacts in the early phase of star formation as they redistribute magnetic flux and suppress angular momentum transport by magnetic fields. We perform threedimensional nestedgrid radiation magnetohydrodynamic (RMHD) simulations including Ohmic dissipation and ambipolar diffusion. Without these effects, magnetic fields transport angular momentum so efficiently that no rotationally supported disk is formed even after the second collapse. Ohmic dissipation works only in a relatively high density region within the first core and suppresses angular momentum transport, enabling formation of a very small rotationally supported disk after the second collapse. With both Ohmic dissipation and ambipolar diffusion, these effects work effectively in almost the entire region within the first core and significant magnetic flux loss occurs. As a result, a rotationally supported disk is formed even before a protostellar core forms. The size of the disk is still small, about 5 AU at the end of the first core phase, but this disk will grow later as gas accretion continues. Thus the nonideal MHD effects can resolve the socalled magnetic braking catastrophe while maintaining the disk size small in the early phase, which is implied from recent interferometric observations.
Subject headings:
stars: formation — ISM: clouds — ISM: jets and outflows — radiative transfer — magnetohydrodynamics1. Introduction
Circumstellar disks are supposed to form as natural byproducts of star formation processes due to large angular momenta in natal molecular cloud cores. Although extensive observational and theoretical studies have been done, theories of disk formation are not complete yet. A disk plays crucial roles in star formation and evolution since its presence affects the evolution of the central star and its feedback. Moreover, circumstellar disks, or protoplanetary disks, are very important not only in the context of star formation but also of planet formation because they are the sites of planet formation. Therefore, it is of crucial importance to construct a consistent model of circumstellar disk formation in the context of star formation.
The key to understand disk formation is angular momentum redistribution during protostellar collapse. Gravitational torque is important when magnetic fields are weak (Bate, 1998; Saigo et al., 2008; Bate, 2010; Tomida et al., 2010a; Bate, 2011), but typically magnetic fields observed in molecular clouds (Crutcher, 2012; Li et al., 2014a, and references therein) are strong enough to remove angular momenta efficiently and overcome the centrifugal barrier (Mouschovias & Paleologou, 1979, 1980; Basu & Mouschovias, 1994, 1995a, 1995b; Tomisaka, 1998, 2000, 2002; Machida et al., 2005b, a; Commerçon et al., 2010; Tomida et al., 2010b; Machida et al., 2011a; Bate et al., 2014). However, this magnetic braking is in fact too effective and strongly suppresses formation of rotationally supported disks and binary/multiples at least in the early phase of star formation and with the ideal MHD approximation, which is called the magnetic braking catastrophe (e.g. Mestel & Spitzer, 1956; Mellon & Li, 2008; Allen et al., 2003; Li et al., 2011, 2014b) or fragmentation crisis (Hennebelle & Teyssier, 2008).
Obviously this problem cannot be real because many protoplanetary disks and extrasolar planets have been observed, and because the fraction of binaries or multiples is known to be high (e.g. Raghavan et al., 2010; Duchêne & Kraus, 2013). Therefore, there must be something missing in those early simulations which were highly idealized and simplified. This problem has been actively discussed recently and many solutions have been proposed. For example, nonideal MHD effects such as Ohmic dissipation and ambipolar diffusion can redistribute magnetic flux and suppress angular momentum transport (Dapp & Basu, 2010; Machida & Matsumoto, 2011; Machida et al., 2011a; Dapp et al., 2012; Machida & Hosokawa, 2013; Tomida et al., 2013). Misalignment between magnetic fields and the rotational axis can reduce the efficiency of magnetic braking and enable disk formation when magnetic fields are not too strong (Joos et al., 2012; Krumholz et al., 2013). Or simply, a circumstellar disk can form later when the surrounding envelope almost disappears and magnetic braking becomes inefficient (Mellon & Li, 2009; Machida et al., 2011b; Machida & Hosokawa, 2013). Turbulence that ubiquitously exists in star forming clouds is another promising mechanism, but there are various interpretations how it circumvents the magnetic braking catastrophe; diffusion or reconnection enhanced by turbulence can redistribute magnetic fields (Lazarian & Vishniac, 1999; SantosLima et al., 2012, 2013), incoherent magnetic fields are less efficient to transport angular momentum (Seifried et al., 2013, 2014), it can naturally induce misalignment between magnetic fields and rotation (Joos et al., 2013), or simply turbulence carries additional angular momentum locally. Given such many solutions proposed, the magnetic braking catastrophe is not so catastrophic as it sounds, but it still remains as quantitative questions; when and how is a disk formed, and how large (massive) is it?
It is observationally well known that star forming clouds are strongly magnetized. Magnetic fields in dense molecular cloud cores are slightly weaker than or comparable to the critical field strength required to support those cloud cores (Crutcher, 2012; Li et al., 2013). Observations of circumstellar disks around young stellar objects also have been actively carried out. Maury et al. (2010) performed interferometric observations toward Class0 sources, and compared them with numerical simulations with and without magnetic fields. They did not detect largelyextended circumstellar disks, suggesting that the effects of magnetic fields are significant in formation and early evolution of circumstellar disks. Tobin et al. (2012, 2013) (see also Sakai et al., 2014) claimed that they found a large circumstellar disk with a radius of around a very young Class0 object L1527 IRS, one of the beststudied examples. While this result was striking, Ohashi et al. (2014) observed the same object with higher resolution and sensitivity using the Atacama Large Millimeter/submillimeter Array (ALMA), and concluded that the disk radius should be about or smaller than . Clearly more detailed and extensive observations are needed, but overall, these observations imply that young circumstellar disks should be small, while they could form in the early phase of star formation and should grow later.
Thus magnetic fields play crucial roles in star formation, and therefore their effects must be carefully clarified and quantified. Because the ionization degree in star forming clouds are extremely low, redistribution of magnetic fields by nonideal MHD effects is naturally expected and has a significant impact on angular momentum transport (Wardle & Koenigl, 1993; Fiedler & Mouschovias, 1993; Basu & Mouschovias, 1994; Ciolek & Mouschovias, 1994; Desch & Mouschovias, 2001; Nakano et al., 2002; Tassis & Mouschovias, 2007; Mellon & Li, 2009; Kunz & Mouschovias, 2009, 2010; Li et al., 2011; Dapp et al., 2012). Nonideal MHD effects are also important in the context of the magnetic flux problem, which means that typical magnetic flux in initial molecular cloud cores is far larger than that in formed stars by orders of magnitude and therefore substantial flux redistribution is required (e.g. Paleologou & Mouschovias, 1983; Braithwaite, 2012). In this work, we investigate the nonideal MHD effects and disk formation during protostellar collapse. For this purpose, we perform threedimensional RMHD simulations of protostellar collapse including nonideal MHD effects. In Tomida et al. (2013) (hereafter Paper I), we found that Ohmic dissipation can redistribute magnetic flux from the dense region and enable formation of rotationally supported disks around protostellar cores. In this paper, we extend our simulations to include ambipolar diffusion, which is more efficient in the lower density region. This paper is organized as follows. The equations and numerical methods are explained in Section 2, and the models used in the simulations are described in Section 3. The results of the numerical simulations are presented in Section 4. Section 5 is devoted to discussion and conclusions are summarized in Section 6. In the Appendix, we describe the newly implemented solver for the nonideal MHD effects and demonstrate its validity.
2. Methods
2.1. Basic Equations
In order to study protostellar collapse, many physical processes must be properly taken into account, such as MHD, selfgravity, radiation transfer, chemistry and nonideal MHD effects. In Paper I, we only considered Ohmic dissipation as nonideal MHD effects. In this work, we extend our threedimensional (3D) nestedgrid simulations to include ambipolar diffusion, which is caused by decoupling between neutral and charged particles. We adopt the single fluid approach assuming strong coupling (Mac Low et al., 1995; Duffin & Pudritz, 2008; Masson et al., 2012). Note that the Hall effect, another possibly important nonideal MHD effect (Li et al., 2011; Braiding & Wardle, 2012), is not included in this work. The governing equations we use are as follows:
(1)  
(2)  
(3)  
(4)  
(5)  
(6)  
(7)  
where are the gas density, gas velocity, magnetic flux density, gas pressure, gas temperature, total gas energy, radiation energy, radiation flux and radiation pressure tensor, respectively. is the speed of light, is the gravitational constant, and is the radiation density constant. From top to bottom, these equations represent conservation of mass, the equation of motion, the induction equation including Ohmic dissipation and ambipolar diffusion, the gas energy equation, the solenoidal constraint, the Poisson’s equation of gravity, and the radiation transfer equations. We adopt the gray flux limited diffusion approximation for radiation transfer (Levermore & Pomraning, 1981; Levermore, 1984). The tabulated Ohmic resistivity and ambipolar diffusion rate are calculated using a chemical network (see 2.3.2). For the Rosseland and Planck mean opacities, and , we combine three tables; Semenov et al. (2003), Seaton et al. (1994) and Ferguson et al. (2005).
2.2. Overview of the Numerical Scheme
In order to achieve the huge dynamic range to resolve star formation processes from the molecular cloud core scale down to the stellar core scale, we use the 3D selfsimilar nestedgrid technique (Yorke & Kaisig, 1995; Ziegler & Yorke, 1997; Machida et al., 2005b, a; Tomida et al., 2013). We solve these equations using the simple operatorsplitting technique, using the same time steps over the whole nestedgrid hierarchy (the socalled shared time stepping). Finer levels are generated at the center of the parent level so that the local Jeans length is resolved at least with 16 cells (Truelove et al., 1997; Commerçon et al., 2008; Joos et al., 2012). We stop the simulations when this condition cannot be satisfied, as unphysical fragmentation may occur. Each level consists of cells, and 14 levels are created by the end of the first core phase. Typical resolution in the first core (35 AU from the center of the cloud, which is covered with level l=12) is or higher. The MHD part is solved using the HLLD approximate Riemann solver (Miyoshi & Kusano, 2005) and the mixed divergence cleaning method (Dedner et al., 2002). The multigrid solver (Matsumoto & Hanawa, 2003) is used for selfgravity. The time step for the MHD and radiation transfer parts is determined by the CourantFriedrichLevy (CFL) condition for the MHD part, and the radiation subsystem is solved implicitly. The implementation of newlyintroduced nonideal MHD effects is described in the Appendix. Detailed descriptions of other parts including microphysical processes are given in Paper I.
Only in the model with ambipolar diffusion, we additionally introduce a density floor where magnetic fields are dominantly strong so that the time step does not become too small. When the gas density is so low that the plasma beta gets below , the gas density is increased so that becomes while the gas temperature and velocity are unchanged. This floor works only in the very low density region in the infalling envelope just above the first core. We monitor the mass added by this density floor and confirm it is very small, about in total, or less than 0.1% of the first core mass. Even with this density floor, the time step is limited by the Alfvén speed in the low density region above the first core in the model with ambipolar diffusion, and as a result an excessive number of time steps are required.
It is important to solve the energy equations including the effects of radiation transfer, although radiation does not drastically affect the evolution in the early phase of star formation compared to the previous simulations using the barotropic approximation (e.g. Bate, 1998; Saigo et al., 2008; Machida et al., 2008a, b). Radiation heating and cooling are only crudely treated in the barotropic approximation, and heating by the nonideal MHD effects is completely neglected. In reality, the thermal evolution varies depending on many factors, but the barotropic approximation simply assumes that the gas temperature is a function of the local gas density. Therefore, it cannot reproduce the realistic evolution and may affect the results artificially. For example, the stability of the first core depends on the gas temperature especially when it is gravitationally unstable (Stamatellos & Whitworth, 2009). Also, the heating by the nonideal MHD effects affects the structure of the first core (see Paper I), and the rates of the nonideal MHD effects are sensitive to the gas temperature. Thus it is crucial to treat both dynamics and thermodynamics consistently.
2.3. Microphysics
Update on the EOS
We update the equation of state from Paper I. Now the partition function for vibration of molecular hydrogen is
(8) 
where is the excitation temperature of vibration of molecular hydrogen
Ohmic Resistivity and Ambipolar Diffusion Rate
The Ohmic resistivity and ambipolar diffusion rate are calculated based on the method described in Nakano et al. (2002) and Okuzumi (2009). This model solves a chemical reaction network of and dust grains, and derive equilibrium solutions. In addition, we consider thermal ionization of potassium (K), which produces a sudden increase of the ionization degree and sharp cutoff in the rates of the nonideal MHD effects around . We adopt the reaction rates from the UMIST RATE 2012 database (McElroy et al., 2013). Dust grains can be charged from 14 to +14. We assume that the dust to gas ratio is 0.01, the dust density is , and the size of the dust particles is . The Ohmic resistivity in this work is higher than that in Paper I in the low density region, by about a factor of 3 in most regions except for the sharp increase around (see Figure 3 below). In Paper I, this rise occurs around a higher density, . These differences result in larger angular momenta remaining in the first cores, but the results are still qualitatively similar.
We simply adopt a constant ionization rate , assuming cosmic rays are the dominant ionization source and neglecting their shielding. In dense regions where cosmic rays cannot penetrate, this assumption overestimates the ionization rate and hence underestimates the diffusion rates. The ionization rate in such regions is generally determined by decay of radioactive elements. At least for our solar system, there is meteoritic evidence that shortlived radionuclides such as existed when the solar nebula formed (Adams, 2010). As discussed in Paper I, these shortlived nuclides yield an ionization rate of (Umebayashi & Nakano, 2009). In this case, neglecting cosmicray shielding would have only a moderate impact. By contrast, the ionization rate in dense regions can be extremely low if only longlived radio active nuclides such as exist (Umebayashi & Nakano, 2009; Kunz & Mouschovias, 2009, 2010; Dapp et al., 2012). Because the lifetime of is short and is comparable to the time scale of star formation, its abundance would vary in each system. Therefore this is an environmental parameter rather than an uncertainty. In this work, we aim to demonstrate that disk formation is possible even with the conservative case with the high ionization rate (i.e. the resulting nonideal MHD effects are weak) as the first step.
Another significant simplification in our model is that dust particles have a single population of , neglecting their distribution and evolution. While dust grains can grow in collapsing clouds, they also can be spattered at shocks. On the other hand, dust grains evaporate at high temperatures. Although the effects of dust evaporation are taken into account in the opacities (Semenov et al., 2003), they are not included in our rates. Specifically, water ices evaporate around , while silicates evaporate around . The latter has no significant impact on the rates because thermal ionization of potassium occurs below that temperature, and magnetic fields and gas are well coupled there. The evaporation of water ices, however, can be important but its effect strongly depends on dust properties such as composition and structures, which are highly uncertain.
Thus the rates of the nonideal MHD effects have large uncertainties, and they would vary in different environments. Therefore, despite the detailed calculation of the rates, our simulations should be considered as an experiment with a typical model, and broad parameter surveys must be performed in the future.
3. Models
Model  OD  AD  RSD Formation  

I  2.4  20  3.8  N  N  770  0.03  2,500  Not Formed  – 
O  2.4  20  3.8  Y  N  1,050  0.04  7,200  After Second Collapse  1 
OA  2.4  20  3.8  Y  Y  2,750  0.075  37,000  Before Second Collapse  5 
Notes. The first five columns are the initial model parameters: the angular velocity, the magnetic field strength, the normalized masstoflux ratio, whether Ohmic dissipation and ambipolar diffusion are introduced. Other parameters are common: , , , and . The last five columns indicate the results of the simulations: the first core lifetime, the first core mass, the first core masstoflux ratio (not normalized), remarks on formation of rotationally supported disks, and the disk radius at the ends of the simulations. See the text for details.
We adopt the same initial conditions as in Paper I. An unstabilized BonnorEbert (hereafter BE, Bonnor, 1956; Ebert, 1955) sphere with m=2 perturbation is used as the initial density profile:
(9) 
where is the the critical BE sphere with and , the radius of the critical BE sphere, the amplitude of the initial density enhancement, the amplitude of regularized perturbation, respectively. We adopt and . The mass and radius of this cloud are and . The initial free fall time at the center of the cloud is .
We initialize the cloud with solidbody rotation and uniform magnetic fields both aligned to axis. The angular rotation speed is or . The magnetic field strength is and the corresponding masstoflux ratio normalized by the critical value of stability for a uniform sphere is where and (Mouschovias & Spitzer, 1976). These parameters are the same to those in the fast rotating models RF and IF in Paper I.
The boundary conditions are also similar to those in Paper I. All the variables except magnetic fields outside the initial BE sphere are fixed to their initial values. Periodic boundary conditions are adopted for magnetic fields in order to avoid the divergence error. Because the gas outside the BE sphere is fixed, the magnetic fields are also anchored to the ambient gas. The boundary condition for selfgravity is set by the multipole expansion (Matsumoto & Hanawa, 2003). These model setups are chosen in order to model isolated dense molecular cloud cores.
We calculate three models; Model I is the ideal MHD approximation, Model O includes Ohmic dissipation only, and Model OA includes both Ohmic dissipation and ambipolar diffusion. The model parameters are summarized in Table 1. While Model I evolves without any trouble until the central temperature exceeds , the other two models are more troublesome. As a result of suppression of angular momentum transport, rotationally supported disks are formed in these two models. These disks become gravitationally unstable after the second collapse in Model O and during the second collapse in Model OA. Because of the limitation of the nestedgrid code, we only can refine the center of the computational domain, and therefore we cannot follow further evolution when fragmentation occurs outside the finest grid. Thus we stop the calculations around in Model O and in OA. For this reason, we focus on the first core phase in this work. Larger grids or more flexible adaptive mesh refinement (AMR) is required for further studies after fragmentation occurs.
4. Results
4.1. Overview of Evolution
Here we summarize the evolutions of the three models. The time evolution of the density and temperature at the center of the clouds as well as the first core masses are shown in Figure 1. Figure 2 is the evolution tracks of the central gas element in the  plane, and Figure 3 shows the rates of Ohmic dissipation and ambipolar diffusion at the center of the cloud in Model OA. The profiles of the density, temperature, velocities and various forces at the end of the first core phase (we simply define it as although the effective adiabatic index gets below the critical value slightly below this temperature) are shown in Figure 4. Model I is essentially the same as Model IF in Paper I, and Model O is similar to RF in Paper I.
Initially the clouds collapse isothermally because radiation cooling is very effective, but the gas starts to evolve almost adiabatically when the central density exceeds about . Then the temperature quickly rises and quasihydrostatic objects, the first cores, are formed (Larson, 1969; Masunaga & Inutsuka, 2000). The formation of the first core is slightly earlier in Model OA, because ambipolar diffusion redistributes the magnetic flux before the gas becomes dense enough to form a first core, and it slightly weaken deceleration of the infall motion by magnetic fields. However, this is not significant because the dynamical collapse is much faster than the redistribution of the magnetic fields, and the difference is only 1% compared to the time from the beginning of the collapse to the formation of the first core. It should be noted that this difference depends on the stability of the initial condition. If the initial condition is sufficiently unstable (super critical), ambipolar diffusion has only minor impact on the initial isothermal collapse (see also Chen & Ostriker, 2014; Heitsch & Hartmann, 2014).
Substantial differences arise during the first core phase. The first core evolves more slowly in the nonideal MHD models, especially with ambipolar diffusion. The nonideal MHD effects redistribute magnetic fluxes from the first cores and suppresses magnetic angular momentum transport. The additional rotational support changes the structures of the first cores and extends their lifetimes. Here we define the first core lifetime simply as the time between the epochs when the central temperature exceeds 100K and 2,000K. The lifetimes of the first cores are about 770 years in Model I, 1,050 years in O, but 2,750 years in OA. The first core in OA is qualitatively different and lives significantly longer than the others. While the accretion rates are essentially the same in all the models (initially and gradually decrease later), the final masses are significantly different, about in Model I, in O, and in OA. The first core in the ideal MHD model is essentially the same with spherical first cores without rotation, while the first core in Model O is slowly rotating but still looks almost spherical as rotational support is not dominant (Figure 4, see also Paper I). The first core in Model OA becomes disklike and is clearly supported by rotation, although contribution of the gas pressure is still substantial. The nonideal MHD effects also provide additional heating as in eq.(4), and the gas element follows a hotter (i.e. higher entropy) evolution track in Model O, and it is more prominent in OA (Figure 2). These variations in the thermodynamic evolutions cannot be reproduced with the barotropic approximation. Because the hotter gas has higher pressure and can support more mass, this also contributes to the larger masses and longer lifetimes in the nonideal MHD models.
When the temperature reaches about , all the dust components evaporate and the opacities drop drastically (Figure 2, see also Paper I). This effect is not obviously visible in the evolution track of Model I because it collapses very quickly. Model O, on the other hand, experiences considerable radiative cooling after dust evaporation and its evolution track merges to that of the ideal MHD model. In Model OA, the dust evaporation affects the thermal evolution more prominently because it evolves more slowly and therefore has more time to cool (see also Tomida, 2014). It evolves along the dust evaporation temperature for a while
Dissociation of molecular hydrogen becomes significant when the temperature at the center of the clouds reaches about . Because of this strongly endothermic reaction, the gas pressure fails to maintain balance with the gravity and starts to collapse again. This second collapse is fairly violent in Model I and O. They dynamically collapses almost at the time scale of freefall. After the second collapse completes, another hydrostatic objects, protostellar cores are formed. While the protostellar core in Model I is almost spherically symmetric, that in Model O is rotating and a disk is simultaneously formed (see also Paper I). Because the resistivity in this work is higher than that in Paper I, the disk is larger but its radius is still small at the end of the simulation, which is only . On the other hand, Model OA evolves again differently; because of the centrifugal force, the second collapse occurs rather gradually. Radiation cooling is still effective even in this phase and the gas evolution follows the colder track, deviating from the adiabatic evolution (Figure 2). As mentioned above, we stop the calculations of Model O and OA earlier because fragmentation occurred.
The evolution of the diffusion rates and corresponding magnetic Reynolds numbers evaluated at the center of the cloud in Model OA is shown in Figure 3. The Ohmic dissipation rate in Model O behaves almost the same. The ambipolar diffusion rate is higher than the Ohmic dissipation rate in almost all the regions, and Ohmic dissipation becomes dominant over ambipolar diffusion only in the dense region where . The Ohmic dissipation rate is very low in the low density region, and it becomes active only in the dense region within the first core (). The ambipolar diffusion rate is much higher than the Ohmic dissipation rate in the low density region, but it does not change the dynamics substantially until the first core is formed. Ambipolar diffusion works effectively in the broader range than the Ohmic dissipation (), removes magnetic flux and suppresses angular momentum transport more significantly. Both the Ohmic dissipation and ambipolar diffusion rates drop sharply when the gas temperature exceeds about (or ) because the gas and magnetic fields recouple by thermal ionization of potassium (K).
4.2. Outflows
First, we compare the outflows driven around the first cores. Vertical cross sections of the density and temperature in this scale are presented in Figure 5. These outflows are driven by the magnetocentrifugal force (Blandford & Payne, 1982) and remove angular momenta from the central objects efficiently (e.g. Tomisaka, 1998, 2002; Tomida et al., 2013). The driving regions are located around , corresponding to the outermost peaks of the rotational velocity in Figure 4. The gas density at these radii is still low and neither Ohmic dissipation nor ambipolar diffusion works strongly yet. The difference of the driving radii between the models is simply originated from the ages of the systems, because accreting angular momenta get larger as the systems evolve.
The sizes of the outflows are significantly different between the models. Before the end of the first core phase, the outflows extend vertically about 70 AU in Model I, 80 AU in O, and 170 AU in OA. However, the structures and velocities of the outflows are similar, about in all the cases. This is because the nonideal MHD effects, both Ohmic dissipation and ambipolar diffusion, only work in the high density region within the first core, and do not affect the driving regions which are located outside the first cores as discussed above (Figure 4). Therefore, the larger sizes in the nonideal MHD models are solely consequences of the longer first core lifetimes by rotational support.
While the outflow structures are quite similar, the temperature distribution in this scale is prominently warmer in Model OA than the other models. This is because the first core in OA is more massive and hotter, and heats up the surrounding material by radiation more strongly. Because of the longer lifetime and higher luminosity, the probability of finding a first core by observations would be higher than previously expected, although it should be still very rare because the lifetime is still much shorter than the timescale of the whole star formation process.
Another noticeable difference in the outflows is that a polar cavity with a very low density is formed in Model OA. The gas near the pole is removed by accretion and the centrifugal force, rather than expelled by the outflow itself. The opening angle of this cavity is still very narrow in the early phase, but it will expand as the gas in the outer envelope with higher angular momenta will accrete later. Although it is not very likely, this cavity might enable us to observe the central object more directly if it is in a faceon configuration.
4.3. First Cores and Second Collapse
As we can see in Figure 3, the nonideal MHD effects are most active in the dense region within the first cores and significant magnetic flux loss occurs where the Reynolds number gets below unity. As a result, much more prominent differences between the models arise in the first core evolution and structures. Figure 6 shows evolution of the masses, angular momenta and masstoflux ratios of the first cores. In order to produce these graphs, we define the first core by a simple criterion; where the gas density is higher than a critical density
Ideal MHD Model
The first core in Model I is essentially the same as that in Model IF of Paper I. Magnetic braking is so efficient that almost all the angular momentum in the first core is removed, and the first core is virtually spherically symmetric and not rotating (Figures 4 and 7). The first core mass and radius at the end of the first core phase are about and 3 AU. The first core and outer pseudodisk exhibit strong warping, which is driven by the magnetic interchange instability (Spruit et al., 1995; Stehle & Spruit, 2001; Krasnopolsky et al., 2012). The angular momentum in the first core decreases as it evolves. The two sharp “wedges” are induced by strong perturbation from the magnetic interchange instability, indicating that the first core angular momentum is very small (essentially zero) and is easily dominated by the external warping. Because the thermal pressure and rotational support are smaller in this model than in the others, and the magnetic pressure is not significant, the first core mass is the smallest in this model.
To see the degree of magnetization, we look into the masstoflux ratio of the first cores. The magnetic flux threading a first core is measured by integrating the vertical magnetic flux density over the midplane as follows:
(10) 
where is the surface area of a discretized cell, and the summation is taken over the midplane. This formula captures the total magnetic flux correctly when the midplane is the plane of symmetry. In Model I, the symmetry breaks down by the magnetic interchange instability after the central density exceeds a few , while the other two models maintain the symmetry during the first core phase. Therefore, the masstoflux ratio in Model I after the sharp increase is not accurate, and it only gives the upper limit. Instead, it is more appropriate to assume this value is almost unchanged beyond this point (the reddotted line in Figure 6) because the time after this is very short and the first core acquires only a little mass (see Figure 1). Taking this point into account, the first core magnetization remains almost constant during the first core phase in the ideal MHD model.
In the large scale, the magnetic fields look like the socalled hourglass shape. The field lines are tangled near and within the first core because of the magnetic interchange instability (Figure 10). Because the first core is not rotating, the field lines are not wound up.
After the second collapse begins, the gas collapses dynamically and a protostellar core is formed very quickly in a few years. The protostellar core is also virtually spherical because almost no angular momentum exists in the protostellar core (at least shortly after its formation). The protostellar core properties are essentially the same as a spherical model without rotation and magnetic fields (Paper I).
NonIdeal MHD Model with Ohmic Dissipation Only
Model O is very similar to Model RF in Paper I, although the Ohmic dissipation rate in this work is higher than the one used in Paper I. Distribution of dissipative regions are shown in Figure 11 (the top row) using the magnetic Reynolds number where is the local Jeans length. Because Ohmic dissipation works only in the high density region within the first core, the first core properties (Figure 6) are identical to those in Model I when it is formed. As the first core evolves, the magnetic flux is redistributed outwardly by Ohmic dissipation. The gas and magnetic fields recouple in the central region after the thermal ionization of potassium begins, but the outer region remains resistive and the masstoflux ratio increases monotonically. At the end of the first core phase, the masstoflux ratio is about three times higher than in Model I. Angular momentum transport by magnetic fields is significantly suppressed, and the first core acquires a large angular momentum, which also monotonically increases. As a result, the first core has considerable rotation and evolves slightly more slowly than Model I (Figure 1). However, the rotation is not fast enough to support the first core and its structure remains almost spherical (Figures 4 and 8). At the end of the first core phase, the mass and radius are about 0.04 and 5 AU, which are larger than the ideal MHD case because of additional support by rotation and heat produced by Ohmic dissipation (Figure 2, see also Paper I). However, these differences in the first cores are rather minor. Thus, roughly speaking, Ohmic dissipation does not have a drastic impact on the first core structure and evolution.
The magnetic fields in the outer region of the first core are straightened out by the strong Ohmic resistivity. On the other hand, the magnetic fields are tightly wound up and amplified by the shearing rotation in the hot recoupled region near the center of the first core (Figures 10 and 11). The magnetic pressure in this region is not strong enough to launch outflows during the first core phase. Later in the protostellar core phase, strong toroidal magnetic fields are produced by the faster rotation in the disk and drive fast jets by the magnetic pressure (see Paper I for the detailed discussion about the protostellar cores).
The second collapse occurs dynamically at the beginning as in Model I because the rotational support is not significant. However, since the gas spins up as it collapses, a rotationally supported disk forms around the protostellar core essentially in parallel. The disk size is very small, only about at the end of the simulation (when ). At the same time, fast bipolar jets are launched by the magnetic pressure of the toroidal magnetic fields amplified by the shearing rotation (Paper I, see also Machida et al., 2008a). This disk is larger than that in Paper I as a consequence of the higher Ohmic resistivity, but qualitatively the same. In this sense, Ohmic dissipation resolves the magnetic braking catastrophe and enables formation of a rotationally supported disk in the early phase of star formation. We stop the simulation because the Jeans condition cannot be satisfied beyond this point due to the limitation of the simulation code.
NonIdeal MHD Model with Ohmic Dissipation and Ambipolar Diffusion
The first core in Model OA with both Ohmic dissipation and ambipolar diffusion is qualitatively different from the other two cases. Ambipolar diffusion is more effective than Ohmic dissipation especially in the low density gas (Figure 3), and suppresses angular momentum transport more significantly. The first core is already less magnetized from the beginning, and it loses the magnetic flux continuously during the first core phase (Figure 6). Figure 11 indicates that the nonideal MHD effects work strongly in almost the whole first core, especially in the early phase (see also Susa et al., 2015). At the end of the first core phase, the masstoflux ratio is higher than the ideal MHD model by more than one order of magnitude.
Angular momentum transport by magnetic fields is strongly suppressed in such a diffusive region, and substantial angular momentum remains within the first core. As a result, the first core becomes a disk supported by the centrifugal force, although the gas pressure is still substantially contributing (Figure 4). The rotational profile is not Keplerian yet (i.e. the rotation profile is not ) because of this gas pressure support and the selfgravity. The disk radius is initially about 6 AU when it is formed, gradually shrinks to , and then remains almost constant, because magnetic braking is still working outside the resistive region and regulates angular momentum accretion. Although the disk remains small until the end of the simulation, it accumulates large angular momentum as the gas accretes from the envelope. Eventually it acquires almost one order of magnitude larger angular momentum than that in Model O, while the mass ( at the end of the first core) is larger only by a factor of two. This fast rotation is sufficient to the trigger the gravitorotational instability (Toomre, 1964; Bate, 1998; Saigo & Tomisaka, 2006; Saigo et al., 2008; Tsukamoto & Machida, 2011), and the disk becomes nonaxisymmetric in the late phase (Figure 9). Angular momentum transport by gravitational torque via nonaxisymmetric structures is usually less effective than magnetic braking, but it becomes important where the magnetic fields are weakened by the nonideal MHD effects. It is naturally expected that this disk will grow later and the rotational support will be more significant as the gas with higher angular momenta will accrete from the outer region.
The magnetic flux expelled by the nonideal MHD effects piles up outside the diffusive region (), i.e. the first core. In this situation, the socalled magnetic wall can be formed, where the magnetic pressure strongly decelerates the infalling gas and possibly form a shock (e.g. Li & McKee, 1996; Tassis & Mouschovias, 2005a, b, 2007; Kunz & Mouschovias, 2010). However, such a structure is not prominent in this model, although the accretion velocity is indeed slower than the freefall velocity due to the magnetic fields (Figure 4). The outermost deceleration (at ) is formed by interacting with the bipolar outflow (see Figure 5), the second one (at ) is a centrifugal barrier, and the innermost one (at ) is the first core shock. This is likely because our rates of the nonideal MHD effects are not high enough, or because the simulation duration is not long enough to accumulate sufficient magnetic flux.
As discussed before, this model suggests that first cores are more readily observed because of the higher luminosity and longer lifetime. Another interesting observational aspect is its variability. While mass and angular momentum transport by magnetic fields are almost steady and smooth, the gravitational torque via nonaxisymmetric structures is unsteady and episodic (see also Vorobyov & Basu, 2006, 2010; Stamatellos et al., 2011). When the disk becomes gravitationally unstable by mass accretion or radiation cooling, nonaxisymmetric structures like spiral arms are formed and transport mass and angular momentum, then the disk is stabilized again unless it is too unstable and fragments. This stochastic accretion causes the first core flickering. Its time scale is corresponding to the orbital time scales, ranging from a few years to a few decades. We will discuss its observational properties in a separate paper in the near future.
Similar to the case with Ohmic dissipation only, the magnetic field lines in the outer region of the first core are strongly straightened out, running almost vertically, by strong Ohmic dissipation and ambipolar diffusion (Figures 10 and 11). The field lines are wound up tightly near the center of the first core where the magnetic fields are coupled with the gas, but again the magnetic pressure is not dominant during the first core phase and the fields behave passively in this phase.
Although we stop the simulation soon after the second collapse begins (when ), we can see that the evolution during the second collapse is already different from the other two models. Because it is supported by the centrifugal force, the gas does not collapse dynamically even when it loses the gas pressure support. It rather collapses gradually with a time scale much longer than the free fall. Unfortunately, we have to stop the simulation at this point because the Jeans condition is violated and fragmentation occurred. We cannot tell but only can speculate how the system will evolve later, but it is likely that a compact binary or multiple system is formed (Machida et al., 2008b).
To summarize, ambipolar diffusion is so effective and drastically changes the first core structure. It removes magnetic flux from the first core, suppresses angular momentum transport, and enables early formation of rotationally supported disks, even before formation of the protostar. The disk is massive enough to become gravitationally unstable and becomes nonaxisymmetric. On the other hand, the disk size remains small in the early phase, contrary to purely hydrodynamic models without magnetic fields which yields very large disks supported by rotation even in the first core phase (Bate, 1998, 2010; Saigo et al., 2008; Tsukamoto & Machida, 2011; Tomida, 2014). Thus, these nonideal MHD effects can solve the magnetic braking catastrophe and maintain consistency with observations of young stellar objects (e.g. Maury et al., 2010).
5. Discussion
5.1. Comparison with the Precedent NonIdeal MHD Simulations
Our results seem to be inconsistent with the earlier nonideal MHD simulations by Li et al. (2011) (see also Mellon & Li, 2009), in which no rotationally supported disk was formed even with ambipolar diffusion. Because there are many differences in the simulation setups including the initial conditions, comparison is not trivial but at least the diffusion rates are similar. We suspect that the most important difference between their simulations and ours is the boundary condition at the center of the cloud. They have a relatively large hole with a radius of 6.7 AU, while we do not have any boundary there but resolve the smallscale structure directly using the nestedgrid technique. In fact, because the disks formed in our simulations are small, it is not inconsistent with their simulations, at least in the early phase of protostellar collapse.
Our calculations are expensive and we cannot follow the longterm evolution, but this high resolution is required to resolve the small first core disk. We believe it is crucial to resolve the first core because the nonideal MHD effects work effectively and a rotationally supported disk is formed there. Otherwise, these nonideal MHD effects have only minor influence on the large scale evolution. Indeed, Machida et al. (2014) demonstrated that the sink particle or inner hole has a significant impact on the results. We propose that the effects of these numerical tricks can be estimated by measuring the angular momentum accreted into the inner boundary or sink particles.
5.2. Perspectives of Future Observations
As discussed in the introduction, recent interferometric observations of Class0 objects, the youngest protostars including L1527 IRS, suggest that circumstellar disks in the early phase of star formation should be small, i.e. , while they can be formed in the early phase (Maury et al., 2010; Ohashi et al., 2014). These observational results do not favor previous hydrodynamic simulations of protostellar collapse without magnetic fields which tend to produce large rotationally supported disks in the early phase (e.g. Bate, 1998, 2010; Saigo et al., 2008; Tsukamoto & Machida, 2011). In contrast, our model with both ambipolar diffusion and Ohmic dissipation is consistent with this picture. In order to test the scenario, more observations of very young stellar objects with higher resolutions that can sufficiently resolve the circumstellar disks are required. In particular, if we could directly observe a first core or at least its remnant shortly after the second collapse, it would be a critical test for our models. If we could detect a compact rotating disk in a collapsing cloud but associated with no significant stellar signature (e.g. near infrared emission or high velocity jets), it would be crucial evidence indicating that a rotationally supported disk is already formed in the first core phase.
From numerical simulations (e.g. Boss & Yorke, 1995; Tomida et al., 2010a; Tomisaka & Tomida, 2011; Saigo & Tomisaka, 2011; Commerçon et al., 2012), a first core is expected to (1) have a spectral energy distribution like a lowtemperature gray body without a hot infrared stellar emission, (2) be a compact, condensed core which may exhibit a signature of infall motion, and (3) be associated with slow compact molecular outflows, but no fast jets. Although first cores should be very rare due to the short lifetime that only one out of a few hundreds of molecular cloud cores are expected, some candidates are already reported (Chen et al., 2010; Enoch et al., 2010; Chen et al., 2012; Pezzuto et al., 2012; Hirano & Liu, 2014; Tokuda et al., 2014). Among them, L1451mm (Pineda et al., 2011) is one of the best examples which satisfies all the above criteria. We expect that direct observation of a first core or its remnant with high resolution interferometers, especially ALMA, will be achieved in the near future.
5.3. NonIdeal MHD Effects
While the nonideal MHD effects can avert the magnetic braking catastrophe, the magnetic flux problem may remain a serious problem. As shown in Figure 6, the nonideal MHD effects increase the masstoflux ratio of the first core but only by one order of magnitude. This is insufficient to explain the difference of magnetization between molecular cloud cores and protostars. This is partly because we do not consider shielding of cosmic rays (see also Dapp et al., 2012; Padovani et al., 2014), while the ionization rate in the dense region would largely vary depending on the abundance of radioactive nuclei as discussed in 2.3.2. Possibly, higher nonideal MHD rates or other mechanisms to accelerate diffusion (e.g. turbulence) are required. If we use higher rates, the disk size will be also larger. Or, we speculate that further redistribution of magnetic fields from the formed protostar will take place in the later phase (e.g. Braithwaite, 2012).
The nonideal MHD effects play an important role in numerical simulations; they actually help convergence. In the ideal MHD approximation, dissipation of magnetic fields is totally determined by resolution. In protostellar collapse simulations, the midplane is a current sheet where the magnetic fields above and below are antiparallel. In ideal MHD, the current sheet can be extremely thin and convergence is not well defined, at least difficult to achieve. For example, Joos et al. (2012, 2013) reported poor convergence of their ideal MHD simulations with misaligned magnetic fields or turbulence; more than 20 cells per the local Jeans length are required for the misaligned cases while 15 cells seem to be sufficient for the aligned cases. With the nonideal MHD effects, on the other hand, the current sheet has a finite thickness, and therefore convergence is now well defined and can be achieved with a finite resolution. However, it still requires high resolution to resolve the dissipation scales sufficiently because the dissipation scale is very small due to the small diffusivity in the low density region.
Although our simulations include as many physical processes as possible, a possibly important effect missing in this work is the Hall effect. The Hall effect is another nonideal MHD effect which can spin up the gas and help disk formation (Krasnopolsky et al., 2011; Li et al., 2011; Braiding & Wardle, 2012). Because the rate of the Hall effect is as high as the ambipolar diffusion (Wardle & Ng, 1999; Nakano et al., 2002; Li et al., 2011), it is likely important in star and disk formation processes. However, it is not easy to implement in highresolution 3D simulations, especially when we resolve the first core directly where the Hall effect works most actively. This is because it requires very small time steps proportional to where is the rate of the Hall effect. The STS acceleration is not effective for the Hall effect because of its hyperbolic nature. Moreover, unlike Ohmic dissipation and ambipolar diffusion, the Hall effect changes the characteristic speeds of the MHD waves (e.g. whistler waves). Therefore, to capture this effect accurately, it is inadequate in principle to use the operatorsplitting and subcycling techniques for the Hall effect, like we did for Ohmic dissipation and ambipolar diffusion. Instead, we have to solve the MHD part and the nonideal MHD part using the same, very small time steps required by the Hall effect, which will be extremely expensive.
Finally, it should be noted that there are significant uncertainties in the rates of the nonideal MHD effects because of large uncertainties in the underlying dust properties as discussed in 2.3.2. The uncertainties also exist in the dust opacities used in the radiation transfer part. Because these uncertainties have direct impacts on formation and evolution of first cores, the results, especially the disk size, might be easily changed by an order of magnitude. Therefore, this work should be considered as a case study based on simple assumptions on the dust properties. A broad parameter survey is required to clarify the effects of the dust properties such as the structure and size distribution. It will be important and interesting to solve evolution of dust grains, including both growth and destruction, consistently in parallel with RMHD simulations. There are many topics on dust left to be explored for the future works.
6. Conclusions
In order to study the effects of Ohmic dissipation and ambipolar diffusion on protostellar collapse and disk formation, we performed 3D nestedgrid RMHD simulations. We simulated only the early phase of star formation before formation of protostellar cores, but we found that significant differences already arose at this point. We revealed that the nonideal MHD effects qualitatively change the structure and evolution of the first cores. Especially, ambipolar diffusion expands the diffusive region and strongly suppresses angular momentum transport, resulting in early formation of rotationally supported disks. The conclusions are summarized as follows:

Ohmic dissipation works in the high density region and suppresses angular momentum transport. As a result, a very small rotationally supported disk is formed after the second collapse, essentially in parallel with formation of the protostellar core.

Ambipolar diffusion extends the diffusive region to the lower density gas, and suppresses angular momentum transport more effectively. The first core is already supported by rotation significantly from the beginning, although the thermal pressure is still important.

The formed disks remain small in the early phase of star formation and the rotation profile is not Keplerian yet, but they will grow later as gas accretion with larger angular momenta continues.

The first core disk with ambipolar diffusion acquires sufficiently large mass and angular momentum to become gravitationally unstable and form nonaxisymmetric structures. As a result, the accretion becomes naturally episodic. Although this is speculative at this point, it is likely that disk fragmentation occurs and a binary/multiple system will be formed later.

Because of additional support by rotation, the first core is more massive, longerlived, and more luminous in the model with ambipolar diffusion. This increases the probability of first core detections, although it is still rare. Also, the first core can exhibit variability by the gravitational instability with the short time scale on the order of years.
These results can reconcile the observational results and theoretical studies. In previous hydrodynamic simulations of protostellar collapse without magnetic fields, large rotationally supported disks are formed in the early phase of star formation, but those results are not favored by the observations of very young stellar objects. In ideal MHD simulations, on the other hand, the magnetic braking is so efficient that a circumstellar disk is difficult to form, which is referred as the magnetic braking catastrophe. Our results suggest that the nonideal MHD effects can resolve the magnetic braking catastrophe in the early phase of star formation even before a protostellar core is formed, and enable formation of a small disk in the early phase. It is expected that the disk will grow later as larger angular momenta accrete.
In previous studies of disk evolution and planet formation, simplified models of isolated, alreadyestablished protoplanetary disks like the minimum mass solar nebula (Hayashi, 1981) have been often used as the initial and boundary conditions. However, our results suggest that a seed of a circum“stellar” disk can be formed in the early phase even before a protostar is formed, it coevolves with the central protostar, and disk fragmentation can happen while the disk and star are still growing via accretion. This means that we need to consider more realistic situations in order to understand disk evolution and planet formation. Ultimately, we need a consistent scenario of star, disk and planet formation – or stellar system formation.
Appendix A A. Numerical Methods for Ohmic Dissipation and Ambipolar Diffusion
a.1. A.1. Conservative Form
In this appendix, we describe how to update the induction and energy equations related to Ohmic dissipation and ambipolar diffusion. We adopt the operatorsplitting approach and solve these effects separately from the other physical processes. The equations solved in this part are as follows:
(A1)  
(A2)  
We can rewrite eq. (A1) into a conservative form:
(A3) 
where , and are the flux tensors of Ohmic dissipation and ambipolar diffusion, respectively. The fluxes of component in direction are:
(A4)  
(A5) 
Note that the dot product () of a vector and a tensor is defined as summation over the first suffix of the tensor; using the Einstein summation convention, which yields a vector. With these fluxes, considering that these flux tensors are antisymmetric and their diagonal components are zero (), eq. (A2) is also rewritten like:
(A6) 
Again note that the dot product () of a tensor and a vector is summation over the second suffix of the tensor; .
a.2. A.2. Correction Terms for
Magnetic fields are defined at the cell centers in our scheme. In order to satisfy the solenoidal constraint (eq.5), we use the mixed cleaning scheme (Dedner et al., 2002) in the MHD part. By introducing additional correction terms for the nonideal MHD part, we can improve the nature of the equations and reduce the divergence error. For Ohmic dissipation, this is achieved using the method proposed by Graves et al. (2008), adding correction terms proportional to in the diagonal components of the Ohmic dissipation fluxes:
(A7) 
For the ambipolar diffusion part, this is done by introducing additional source terms (Duffin & Pudritz, 2008):
(A8)  
(A9) 
a.3. A.3. Discretization
Discretization of eqs.(A1) and (A2) is straightforward:
(A10)  
(A11)  
where and are the spatial indices of the cell centers in the and direction, and means the location of the cell surface. The Einstein notation is used in eq.(A11); , etc.. The fluxes at the cell surface are calculated using the following relations and their permutations:
(A12)  
(A13)  
(A14)  
(A15)  
(A16) 
For the source terms in eqs.(A8) and (A9), the derivatives at the cell centers are discretized using the following formula and its permutations:
(A17) 
a.4. A.4. NestedGrid Boundaries
At the boundaries between different nestedgrid levels, we have to connect solutions in different levels. Because our scheme uses the conservative form, this is straightforwardly implemented using the following procedure.

Project the solution at from the coarser grids to the finer grid boundaries using linear interpolation.

Calculate fluxes in all the grid levels.

Recalculate the fluxes at the level boundaries in the coarser grid using the sum of the fluxes in the finer grid; e.g. where and are the corresponding cell indices in the coarser and finer grid levels, is the index of the nestedgrid level (larger is finer), and and are the surface areas of the coarser and finer cells, respectively.

Update the solution to in all the grid levels using the corrected fluxes.
Note that no timeinterpolation is required because we use shared time stepping.
a.5. A.5. Time Integration: Super Time Stepping
As in Paper I, we adopt the Super Time Stepping method (STS, Alexiades et al., 1996; Choi et al., 2009; Commerçon et al., 2011) in order to accelerate the calculation when the time step for the nonideal MHD part is much smaller than the MHD part. There are two parameters in STS; which is a small positive parameter controlling the stability and efficiency of the scheme, and which is the number of the sub steps used in the scheme. We use and as in Paper I. Acceleration by a factor of is achieved with these parameters compared to the firstorder explicit scheme.
We use the same time step in all the grid levels. The time step for the nonideal MHD part is estimated from the standard CFL condition for the diffusionlike equations:
(A18) 
where 0.1 is a safety factor. When this time step is significantly smaller than that for the MHD part (), we update the nonideal MHD part using STS and subcycling. If this condition is not satisfied but the time step is still smaller than that for the MHD part (), the standard forward Euler integrator with subcycling is used. If it is larger than the MHD time step (), we update the nonideal MHD part using the time step for the MHD part.
a.6. A.6. Tests: Oblique Cshock Problems
Here we demonstrate the validity of our nonideal MHD solver based on the nonisothermal oblique Cshock tests (Wardle, 1991; Mac Low et al., 1995; Duffin & Pudritz, 2008; Masson et al., 2012). These are simple onedimensional shocktube tests with the nonideal MHD effects.
Although these problems are onedimensional, we make them twodimensional tests by adopting the skewperiodic boundary condition. We rotate the shock tube tests by an angle of . When is a rational number, we can construct the skewperiodic boundary condition which does not introduce any artificial effect. In the direction, we adopt the fixed boundary condition, which may violate the solenoidal constraint but the error is cleaned by the mixed cleaning scheme (Dedner et al., 2002) and does not affect the calculations. In the direction, we copy all the variables from the other side of the computational domain with a shift;
(A19) 
where is the number of cells in the direction. If we set so that is an integer, this boundary condition produces no artificial effect and can be used to run onedimensional shock tube tests on a twodimensional grid.
We compare the results between three configurations: fiducial highresolution aligned models ( and ), highresolution tilted models (, and ), and lowresolution tilted models ( and ). The rotation angle in the tilted models is . Each cell is square shaped; . The resolution in our highresolution aligned model is corresponding to the highest AMR level in Masson et al. (2012).
We adopt the same parameter sets as Masson et al. (2012). For Ohmic dissipation Cshock, the initial left state is , the right is , and the uniform resistivity is used. For ambipolar diffusion, the left is , the right is , and the ambipolar diffusion rate is density dependent; . We run the simulations until steady solutions are obtained. In both calculations, the adiabatic index is , the CFL number is , and the STS time integrator is selected. The results are shown in Figure 12. Overall, we successfully obtained good numerical solutions consistent with the semianalytic solutions. In particular, we found no problem in the multidimensional (i.e. tilted) models.
Footnotes
 affiliation: Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA; tomida@astro.princeton.edu
 affiliation: Department of Physics, The University of Tokyo, Tokyo, 1130033, Japan
 affiliation: Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Meguroku, Tokyo, 1528551, Japan; okuzumi@geo.titech.ac.jp
 affiliation: Department of Earth and Planetary Sciences, Faculty of Sciences, Kyushu University, Hakozaki, Higashiku, Fukuoka 8128581, Japan; machida.masahiro.018@m.kyushuu.ac.jp
 For example, NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101, Release 16a, http://cccbdb.nist.gov/
 Although this behavior is interesting, this might be artificial because we use the tabulated dust opacities as functions of the density and temperature, assuming there is no hysteresis. Specifically, the evaporated dust grains immediately form again once the gas cools down below the evaporation temperature. In reality, it would take some time to reform dust grains and their structure and/or size distribution would be different from those before the evaporation. While this reformation of dust grains is highly uncertain, fortunately, it would not affect the evolution significantly. The gas evolves along the dust evaporation line if the opacities after reformation are higher. On the other hand, it evolves almost isothermally until the gas opacities and density become large enough if the opacities are lower, but the evolution is already close enough to isothermal.
 This value is one order of magnitude higher than the critical density used in Paper I. We have changed the threshold so as to exclude the high density region in the pseudodisk (a flattened disklike structure but not supported by gas pressure or rotation) especially in Model OA. This is why the angular momentum of Model I shown in Figure 6 behaves differently from that in Paper I, but in fact they are consistent if the same threshold is adopted.
References
 Adams, F. C. 2010, ARA&A, 48, 47
 Alexiades, V., Amiez, G., & Gremaud, P. 1996, Communications in numerical methods in engineering, 12, 31
 Allen, A., Li, Z.Y., & Shu, F. H. 2003, ApJ, 599, 363
 Basu, S., & Mouschovias, T. C. 1994, ApJ, 432, 720
 —. 1995a, ApJ, 452, 386
 —. 1995b, ApJ, 453, 271
 Bate, M. R. 1998, ApJ, 508, L95
 —. 2010, MNRAS, 404, L79
 —. 2011, MNRAS, 417, 2036
 Bate, M. R., Tricco, T. S., & Price, D. J. 2014, MNRAS, 437, 77
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
 Bonnor, W. B. 1956, MNRAS, 116, 351
 Boss, A. P., & Yorke, H. W. 1995, ApJ, 439, L55
 Braiding, C. R., & Wardle, M. 2012, MNRAS, 422, 261
 Braithwaite, J. 2012, MNRAS, 422, 619
 Chen, C.Y., & Ostriker, E. C. 2014, ApJ, 785, 69
 Chen, X., Arce, H. G., Dunham, M. M., et al. 2012, ApJ, 751, 89
 Chen, X., Arce, H. G., Zhang, Q., et al. 2010, ApJ, 715, 1344
 Choi, E., Kim, J., & Wiita, P. J. 2009, ApJS, 181, 413
 Ciolek, G. E., & Mouschovias, T. C. 1994, ApJ, 425, 142
 Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371
 —. 2010, A&A, 510, L3
 Commerçon, B., Launhardt, R., Dullemond, C., & Henning, T. 2012, A&A, 545, A98
 Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35
 Crutcher, R. M. 2012, ARA&A, 50, 29
 Dapp, W. B., & Basu, S. 2010, A&A, 521, L56
 Dapp, W. B., Basu, S., & Kunz, M. W. 2012, A&A, 541, A35
 Dedner, A., Kemm, F., Kröner, D., et al. 2002, Journal of Computational Physics, 175, 645
 Desch, S. J., & Mouschovias, T. C. 2001, ApJ, 550, 314
 Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269
 Duffin, D. F., & Pudritz, R. E. 2008, MNRAS, 391, 1659
 Ebert, R. 1955, Zeitschrift fur Astrophysik, 36, 222
 Enoch, M. L., Lee, J.E., Harvey, P., Dunham, M. M., & Schnee, S. 2010, ApJ, 722, L33
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585
 Fiedler, R. A., & Mouschovias, T. C. 1993, ApJ, 415, 680
 Graves, D. T., Trebotich, D., Miller, G. H., & Colella, P. 2008, Journal of Computational Physics, 227, 4797
 Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
 Heitsch, F., & Hartmann, L. 2014, MNRAS, 443, 230
 Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25
 Hirano, N., & Liu, F.c. 2014, ApJ, 789, 50
 Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128
 Joos, M., Hennebelle, P., Ciardi, A., & Fromang, S. 2013, A&A, 554, A17
 Krasnopolsky, R., Li, Z.Y., & Shang, H. 2011, ApJ, 733, 54
 Krasnopolsky, R., Li, Z.Y., Shang, H., & Zhao, B. 2012, ApJ, 757, 77
 Krumholz, M. R., Crutcher, R. M., & Hull, C. L. H. 2013, ApJ, 767, L11
 Kunz, M. W., & Mouschovias, T. C. 2009, ApJ, 693, 1895
 —. 2010, MNRAS, 408, 322
 Larson, R. B. 1969, MNRAS, 145, 271
 Lazarian, A., & Vishniac, E. T. 1999, ApJ, 517, 700
 Levermore, C. D. 1984, Journal of Quantitative Spectroscopy and Radiative Transfer, 31, 149
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321
 Li, H.b., Fang, M., Henning, T., & Kainulainen, J. 2013, MNRAS, 436, 3707
 Li, H.b., Goodman, A., Sridharan, T. K., et al. 2014a, ArXiv eprints, arXiv:1404.2024
 Li, Z.Y., Banerjee, R., Pudritz, R. E., et al. 2014b, ArXiv eprints, arXiv:1401.2219
 Li, Z.Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180
 Li, Z.Y., & McKee, C. F. 1996, ApJ, 464, 373
 Mac Low, M.M., Norman, M. L., Konigl, A., & Wardle, M. 1995, ApJ, 442, 726
 Machida, M. N., & Hosokawa, T. 2013, MNRAS, 431, 1719
 Machida, M. N., Inutsuka, S.i., & Matsumoto, T. 2008a, ApJ, 676, 1088
 Machida, M. N., Inutsuka, S.I., & Matsumoto, T. 2011a, PASJ, 63, 555
 Machida, M. N., Inutsuka, S.i., & Matsumoto, T. 2011b, ApJ, 729, 42
 —. 2014, MNRAS, 438, 2278
 Machida, M. N., & Matsumoto, T. 2011, MNRAS, 413, 2767
 Machida, M. N., Matsumoto, T., Hanawa, T., & Tomisaka, K. 2005a, MNRAS, 362, 382
 Machida, M. N., Matsumoto, T., Tomisaka, K., & Hanawa, T. 2005b, MNRAS, 362, 369
 Machida, M. N., Tomisaka, K., Matsumoto, T., & Inutsuka, S.i. 2008b, ApJ, 677, 327
 Masson, J., Teyssier, R., MuletMarquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24
 Masunaga, H., & Inutsuka, S.i. 2000, ApJ, 531, 350
 Matsumoto, T., & Hanawa, T. 2003, ApJ, 583, 296
 Maury, A. J., André, P., Hennebelle, P., et al. 2010, A&A, 512, A40
 McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36
 Mellon, R. R., & Li, Z.Y. 2008, ApJ, 681, 1356
 —. 2009, ApJ, 698, 922
 Mestel, L., & Spitzer, Jr., L. 1956, MNRAS, 116, 503
 Miyoshi, T., & Kusano, K. 2005, Journal of Computational Physics, 208, 315
 Mouschovias, T. C., & Paleologou, E. V. 1979, ApJ, 230, 204
 —. 1980, ApJ, 237, 877
 Mouschovias, T. C., & Spitzer, Jr., L. 1976, ApJ, 210, 326
 Nakano, T., Nishi, R., & Umebayashi, T. 2002, ApJ, 573, 199
 Ohashi, N., Saigo, K., Aso, Y., et al. 2014, ArXiv eprints, arXiv:1410.0172
 Okuzumi, S. 2009, ApJ, 698, 1122
 Padovani, M., Galli, D., Hennebelle, P., Commerçon, B., & Joos, M. 2014, A&A, 571, A33
 Paleologou, E. V., & Mouschovias, T. C. 1983, ApJ, 275, 838
 Pezzuto, S., Elia, D., Schisano, E., et al. 2012, A&A, 547, A54
 Pineda, J. E., Arce, H. G., Schnee, S., et al. 2011, ApJ, 743, 201
 Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1
 Saigo, K., & Tomisaka, K. 2006, ApJ, 645, 381
 —. 2011, ApJ, 728, 78
 Saigo, K., Tomisaka, K., & Matsumoto, T. 2008, ApJ, 674, 997
 Sakai, N., Sakai, T., Hirota, T., et al. 2014, Nature, 507, 78
 SantosLima, R., de Gouveia Dal Pino, E. M., & Lazarian, A. 2012, ApJ, 747, 21
 —. 2013, MNRAS, 429, 3371
 Seaton, M. J., Yan, Y., Mihalas, D., & Pradhan, A. K. 1994, MNRAS, 266, 805
 Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2013, MNRAS, 432, 3320
 —. 2014, ArXiv eprints, arXiv:1408.2989
 Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611
 Spruit, H. C., Stehle, R., & Papaloizou, J. C. B. 1995, MNRAS, 275, 1223
 Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 400, 1563
 Stamatellos, D., Whitworth, A. P., & Hubber, D. A. 2011, ApJ, 730, 32
 Stehle, R., & Spruit, H. C. 2001, MNRAS, 323, 587
 Susa, H., Doi, K., & Omukai, K. 2015, ArXiv eprints, arXiv:1501.00087
 Tassis, K., & Mouschovias, T. C. 2005a, ApJ, 618, 769
 —. 2005b, ApJ, 618, 783
 —. 2007, ApJ, 660, 388
 Tobin, J. J., Hartmann, L., Chiang, H.F., et al. 2012, Nature, 492, 83
 —. 2013, ApJ, 771, 48
 Tokuda, K., Onishi, T., Saigo, K., et al. 2014, ApJ, 789, L4
 Tomida, K. 2014, ApJ, 786, 98
 Tomida, K., Machida, M. N., Saigo, K., Tomisaka, K., & Matsumoto, T. 2010a, ApJ, 725, L239
 Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2013, ApJ, 763, 6
 —. 2010b, ApJ, 714, L58
 Tomisaka, K. 1998, ApJ, 502, L163
 —. 2000, ApJ, 528, L41
 —. 2002, ApJ, 575, 306
 Tomisaka, K., & Tomida, K. 2011, PASJ, 63, 1151
 Toomre, A. 1964, ApJ, 139, 1217
 Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179
 Tsukamoto, Y., & Machida, M. N. 2011, MNRAS, 416, 591
 Umebayashi, T., & Nakano, T. 2009, ApJ, 690, 69
 Vorobyov, E. I., & Basu, S. 2006, ApJ, 650, 956
 —. 2010, ApJ, 719, 1896
 Wardle, M. 1991, MNRAS, 251, 119
 Wardle, M., & Koenigl, A. 1993, ApJ, 410, 218
 Wardle, M., & Ng, C. 1999, MNRAS, 303, 239
 Yorke, H. W., & Kaisig, M. 1995, Computer Physics Communications, 89, 29
 Ziegler, U., & Yorke, H. W. 1997, Computer Physics Communications, 101, 54