Non-Ideal RMHD Simulations of Protostellar Collapse

Radiation Magnetohydrodynamic Simulations of Protostellar Collapse: Non-Ideal Magnetohydrodynamic Effects and Early Formation of Circumstellar Disks


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, non-ideal 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 three-dimensional nested-grid 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 non-ideal MHD effects can resolve the so-called 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 — magnetohydrodynamics

1. Introduction

Circumstellar disks are supposed to form as natural by-products 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 extra-solar 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, non-ideal 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). Mis-alignment 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; Santos-Lima et al., 2012, 2013), incoherent magnetic fields are less efficient to transport angular momentum (Seifried et al., 2013, 2014), it can naturally induce mis-alignment 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 Class-0 sources, and compared them with numerical simulations with and without magnetic fields. They did not detect largely-extended 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 Class-0 object L1527 IRS, one of the best-studied 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 non-ideal 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). Non-ideal 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 non-ideal MHD effects and disk formation during protostellar collapse. For this purpose, we perform three-dimensional RMHD simulations of protostellar collapse including non-ideal 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 non-ideal 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, self-gravity, radiation transfer, chemistry and non-ideal MHD effects. In Paper I, we only considered Ohmic dissipation as non-ideal MHD effects. In this work, we extend our three-dimensional (3D) nested-grid 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 non-ideal MHD effect (Li et al., 2011; Braiding & Wardle, 2012), is not included in this work. The governing equations we use are as follows:


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 self-similar nested-grid technique (Yorke & Kaisig, 1995; Ziegler & Yorke, 1997; Machida et al., 2005b, a; Tomida et al., 2013). We solve these equations using the simple operator-splitting technique, using the same time steps over the whole nested-grid hierarchy (the so-called 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 (3-5 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 self-gravity. The time step for the MHD and radiation transfer parts is determined by the Courant-Friedrich-Levy (CFL) condition for the MHD part, and the radiation subsystem is solved implicitly. The implementation of newly-introduced non-ideal 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 non-ideal 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 non-ideal MHD effects affects the structure of the first core (see Paper I), and the rates of the non-ideal 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


where is the excitation temperature of vibration of molecular hydrogen5. This formula is better than the previous one used in Paper I because it does not diverge in the low temperature limit. However, this change has very minor impacts on the results of hydrodynamic simulations because contribution from the vibrational degrees of freedom becomes significant only in a high temperature, and most of hydrogen molecules are already dissociated there. Other parameters are unchanged including the ortho:para ratio of molecular hydrogen, which is 3:1.

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 cut-off in the rates of the non-ideal 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 short-lived radionuclides such as existed when the solar nebula formed (Adams, 2010). As discussed in Paper I, these short-lived nuclides yield an ionization rate of (Umebayashi & Nakano, 2009). In this case, neglecting cosmic-ray shielding would have only a moderate impact. By contrast, the ionization rate in dense regions can be extremely low if only long-lived 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 non-ideal 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 non-ideal 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 mass-to-flux 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 mass-to-flux 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.

Table 1Summary of the initial model parameters and results

We adopt the same initial conditions as in Paper I. An unstabilized Bonnor-Ebert (hereafter BE, Bonnor, 1956; Ebert, 1955) sphere with m=2 perturbation is used as the initial density profile:


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 solid-body rotation and uniform magnetic fields both aligned to -axis. The angular rotation speed is or . The magnetic field strength is and the corresponding mass-to-flux 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 self-gravity is set by the multi-pole 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 nested-grid 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

Figure 1.— Time evolution of the gas density (top) and temperature (middle) at the center of the clouds, and the first core masses (bottom). The temperatures where the dust evaporation completes and where the second collapse begins are also shown with the gray dotted and dash-dotted lines. Note that these temperatures depend on the gas density, but the dependency is weak and the densities where these events happen in these models are close enough (see Figure 2).
Figure 2.— Evolution tracks of the gas elements at the center of the clouds in the - plane. The small inset shows the region near the onset of the second collapse. The gray-dotted line indicates the temperature at which dust grains evaporate completely (Semenov et al., 2003). The evaporation begins below this temperature and the opacity is already reduced. Also, the dash-dotted line roughly shows the temperature where the second collapse begins (i.e. , see also Tomida (2014)).
Figure 3.— The rates and corresponding Reynolds numbers of Ohmic dissipation and ambipolar diffusion at the center of the cloud in Model OA as functions of the gas density. The magnetic Reynolds numbers are calculated assuming that typical velocity and length scale are the free-fall velocity and Jeans length; (see Paper I, Appendix C).
Figure 4.— Radial profiles of the gas density, temperature along the - (in the disk mid-plane, red) and -axes (along the rotational axis, blue), the infall (red) and rotation (blue) velocities along the -axis (in the mid-plane), and the radial accelerations by the various forces (The gas pressure gradient (red), the magnetic pressure gradient (green), the magnetic tension (blue), the centrifugal force (magenta), and the gravity (black). The radiation force is not plotted because it is negligibly small here) in the mid-plane at the end of the first core phase, from top to bottom. From left to right, the different columns are for Model I, O and OA, respectively.

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 quasi-hydrostatic 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 non-ideal MHD models, especially with ambipolar diffusion. The non-ideal 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 disk-like and is clearly supported by rotation, although contribution of the gas pressure is still substantial. The non-ideal 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 non-ideal 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 while6, and the gas gets colder than those in Model I and O.

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 free-fall. 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

Figure 5.— Vertical cross sections of the gas density (top) and temperature (bottom) in the outflow scales ( or for Model I, and or for the rest) at the end of the first core phase (). Projected velocity vectors are overplotted with white arrows. Because the outflows are roughly axisymmetric in the large scale, only vertical cross sections are shown here.

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 magneto-centrifugal 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 outer-most 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 non-ideal 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 non-ideal 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 face-on configuration.

4.3. First Cores and Second Collapse

Figure 6.— Evolution of the masses, angular momenta, and mass-to-flux ratios of the first cores from top to bottom as functions of the gas density at the center of the cloud. The red dotted line in the bottom panel is extrapolation of the mass-to-flux ratio before the magnetic interchange instability develops (see the text). Note that this mass-to-flux ratio is presented in the cgs unit and not normalized.

As we can see in Figure 3, the non-ideal 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 mass-to-flux 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 density7 . Also, the density and temperature cross sections of each model are shown in Figures 7-9. Three epochs are shown in these Figures; just before thermal ionization of potassium (, left), beginning of dust evaporation (, middle) and the end of the first core phase, just after the onset of the second collapse (, right). The magnetic field lines at the end of the first core phase are presented in Figure 10. The properties of the first cores and disks are summarized in Table 1.

Ideal MHD Model

Figure 7.— Vertical (top two rows) and horizontal (bottom two rows) cross sections of the gas density (the first and third rows) and temperature (the second and fourth rows) in the first core scale (l=11 or ) of Model I, when the temperature at the center of the cloud is 600K (left column), 1,400K (middle) and 2,000K (right). is the age of the first core after its formation (when ).

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 pseudo-disk 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 mass-to-flux ratio of the first cores. The magnetic flux threading a first core is measured by integrating the vertical magnetic flux density over the mid-plane as follows:


where is the surface area of a discretized cell, and the summation is taken over the mid-plane. This formula captures the total magnetic flux correctly when the mid-plane 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 mass-to-flux 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 red-dotted 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 so-called 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).

Non-Ideal MHD Model with Ohmic Dissipation Only

Figure 8.— Same as Figure 7 but of Model O.

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 mass-to-flux ratio increases monotonically. At the end of the first core phase, the mass-to-flux 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.

Non-Ideal MHD Model with Ohmic Dissipation and Ambipolar Diffusion

Figure 9.— Same as Figure 7 but of Model OA.

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 non-ideal 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 mass-to-flux 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 self-gravity. 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 gravito-rotational instability (Toomre, 1964; Bate, 1998; Saigo & Tomisaka, 2006; Saigo et al., 2008; Tsukamoto & Machida, 2011), and the disk becomes non-axisymmetric in the late phase (Figure 9). Angular momentum transport by gravitational torque via non-axisymmetric structures is usually less effective than magnetic braking, but it becomes important where the magnetic fields are weakened by the non-ideal 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 non-ideal MHD effects piles up outside the diffusive region (), i.e. the first core. In this situation, the so-called 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 free-fall 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 non-ideal 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 non-axisymmetric 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, non-axisymmetric 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 non-axisymmetric. 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 non-ideal MHD effects can solve the magnetic braking catastrophe and maintain consistency with observations of young stellar objects (e.g. Maury et al., 2010).

Figure 10.— Magnetic field lines (yellow lines) threading the first cores and surrounding pseudo-disks (l=11 or ) at the end of the first core phase. Vertical and horizontal density cross sections are also presented. Note that the density of the field lines does not reflect the strength of the magnetic fields.
Figure 11.— Vertical cross sections of the magnetic Reynolds numbers in the first cores and surrounding envelopes (l=10 or ), when the temperature at the center of the cloud is 600K (left column), 1,400K (middle) and 2,000K (right). We use the gas velocity and local Jeans length to evaluate the Reynolds numbers. The top row shows the Reynolds number calculated from the Ohmic resistivity in Model O. The bottom two rows are Model OA, while the middle row indicates the Ohmic Reynolds number and the bottom row is calculated using the ambipolar diffusion rate. The non-ideal MHD effects work effectively in the blue regions while the red regions are close to ideal.

5. Discussion

5.1. Comparison with the Precedent Non-Ideal MHD Simulations

Our results seem to be inconsistent with the earlier non-ideal 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 small-scale structure directly using the nested-grid 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 long-term 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 non-ideal MHD effects work effectively and a rotationally supported disk is formed there. Otherwise, these non-ideal 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 Class-0 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 low-temperature 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, L1451-mm (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. Non-Ideal MHD Effects

While the non-ideal MHD effects can avert the magnetic braking catastrophe, the magnetic flux problem may remain a serious problem. As shown in Figure 6, the non-ideal MHD effects increase the mass-to-flux 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 non-ideal 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 non-ideal 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 mid-plane is a current sheet where the magnetic fields above and below are anti-parallel. 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 mis-aligned magnetic fields or turbulence; more than 20 cells per the local Jeans length are required for the mis-aligned cases while 15 cells seem to be sufficient for the aligned cases. With the non-ideal 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 non-ideal 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 high-resolution 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 operator-splitting and sub-cycling techniques for the Hall effect, like we did for Ohmic dissipation and ambipolar diffusion. Instead, we have to solve the MHD part and the non-ideal 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 non-ideal 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 nested-grid 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 non-ideal 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:

  1. 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.

  2. 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.

  3. 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.

  4. The first core disk with ambipolar diffusion acquires sufficiently large mass and angular momentum to become gravitationally unstable and form non-axisymmetric 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.

  5. Because of additional support by rotation, the first core is more massive, longer-lived, 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 non-ideal 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, already-established 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 co-evolves 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.

We are grateful to Professor Eve C. Ostriker, Professor James M. Stone, Professor Tomoaki Matsumoto, Dr. Matthew Kunz, Dr. Yasunori Hori and Dr. Yusuke Tsukamoto for fruitful discussions and encouraging comments. We also thank the anonymous referee for the useful comments to improve this paper. This work is partly supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Grants-in-Aid for Scientific Research 23103005, 25887023, 26400224 (SO), 25400232 and 26103707 (MNM). KT is supported by Japan Society for the Promotion of Science (JSPS) Research Fellowship for Young Scientists (). This research used computational resources of the High Performance Computing Infrastructure (HPCI) system provided by the Cybermedia Center at Osaka University and the Cyberscience Center at Tohoku University through the HPCI System Research Project (Project ID: hp140065).

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 operator-splitting approach and solve these effects separately from the other physical processes. The equations solved in this part are as follows:


We can rewrite eq. (A1) into a conservative form:


where , and are the flux tensors of Ohmic dissipation and ambipolar diffusion, respectively. The fluxes of component in -direction are:


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 anti-symmetric and their diagonal components are zero (), eq. (A2) is also rewritten like:


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 non-ideal 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:


For the ambipolar diffusion part, this is done by introducing additional source terms (Duffin & Pudritz, 2008):


a.3. A.3. Discretization

Discretization of eqs.(A1) and (A2) is straightforward:


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:


For the source terms in eqs.(A8) and (A9), the derivatives at the cell centers are discretized using the following formula and its permutations:


a.4. A.4. Nested-Grid Boundaries

At the boundaries between different nested-grid levels, we have to connect solutions in different levels. Because our scheme uses the conservative form, this is straightforwardly implemented using the following procedure.

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

  2. Calculate fluxes in all the grid levels.

  3. 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 nested-grid level (larger is finer), and and are the surface areas of the coarser and finer cells, respectively.

  4. Update the solution to in all the grid levels using the corrected fluxes.

Note that no time-interpolation 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 non-ideal 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 first-order explicit scheme.

We use the same time step in all the grid levels. The time step for the non-ideal MHD part is estimated from the standard CFL condition for the diffusion-like equations:


where 0.1 is a safety factor. When this time step is significantly smaller than that for the MHD part (), we update the non-ideal MHD part using STS and sub-cycling. 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 sub-cycling is used. If it is larger than the MHD time step (), we update the non-ideal MHD part using the time step for the MHD part.

a.6. A.6. Tests: Oblique C-shock Problems

Figure 12.— The results of the oblique C-shock tests with Ohmic dissipation (left) and ambipolar diffusion (right).

Here we demonstrate the validity of our non-ideal MHD solver based on the non-isothermal oblique C-shock tests (Wardle, 1991; Mac Low et al., 1995; Duffin & Pudritz, 2008; Masson et al., 2012). These are simple one-dimensional shock-tube tests with the non-ideal MHD effects.

Although these problems are one-dimensional, we make them two-dimensional tests by adopting the skew-periodic boundary condition. We rotate the shock tube tests by an angle of . When is a rational number, we can construct the skew-periodic 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;


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 one-dimensional shock tube tests on a two-dimensional grid.

We compare the results between three configurations: fiducial high-resolution aligned models ( and ), high-resolution tilted models (, and ), and low-resolution tilted models ( and ). The rotation angle in the tilted models is . Each cell is square shaped; . The resolution in our high-resolution 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 C-shock, 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 semi-analytic solutions. In particular, we found no problem in the multi-dimensional (i.e. tilted) models.


  1. affiliation: Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA;
  2. affiliation: Department of Physics, The University of Tokyo, Tokyo, 113-0033, Japan
  3. affiliation: Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Meguro-ku, Tokyo, 152-8551, Japan;
  4. affiliation: Department of Earth and Planetary Sciences, Faculty of Sciences, Kyushu University, Hakozaki, Higashi-ku, Fukuoka 812-8581, Japan;
  5. For example, NIST Computational Chemistry Comparison and Benchmark Database, NIST Standard Reference Database Number 101, Release 16a,
  6. 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.
  7. 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 pseudo-disk (a flattened disk-like 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.


  1. Adams, F. C. 2010, ARA&A, 48, 47
  2. Alexiades, V., Amiez, G., & Gremaud, P. 1996, Communications in numerical methods in engineering, 12, 31
  3. Allen, A., Li, Z.-Y., & Shu, F. H. 2003, ApJ, 599, 363
  4. Basu, S., & Mouschovias, T. C. 1994, ApJ, 432, 720
  5. —. 1995a, ApJ, 452, 386
  6. —. 1995b, ApJ, 453, 271
  7. Bate, M. R. 1998, ApJ, 508, L95
  8. —. 2010, MNRAS, 404, L79
  9. —. 2011, MNRAS, 417, 2036
  10. Bate, M. R., Tricco, T. S., & Price, D. J. 2014, MNRAS, 437, 77
  11. Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
  12. Bonnor, W. B. 1956, MNRAS, 116, 351
  13. Boss, A. P., & Yorke, H. W. 1995, ApJ, 439, L55
  14. Braiding, C. R., & Wardle, M. 2012, MNRAS, 422, 261
  15. Braithwaite, J. 2012, MNRAS, 422, 619
  16. Chen, C.-Y., & Ostriker, E. C. 2014, ApJ, 785, 69
  17. Chen, X., Arce, H. G., Dunham, M. M., et al. 2012, ApJ, 751, 89
  18. Chen, X., Arce, H. G., Zhang, Q., et al. 2010, ApJ, 715, 1344
  19. Choi, E., Kim, J., & Wiita, P. J. 2009, ApJS, 181, 413
  20. Ciolek, G. E., & Mouschovias, T. C. 1994, ApJ, 425, 142
  21. Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371
  22. —. 2010, A&A, 510, L3
  23. Commerçon, B., Launhardt, R., Dullemond, C., & Henning, T. 2012, A&A, 545, A98
  24. Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35
  25. Crutcher, R. M. 2012, ARA&A, 50, 29
  26. Dapp, W. B., & Basu, S. 2010, A&A, 521, L56
  27. Dapp, W. B., Basu, S., & Kunz, M. W. 2012, A&A, 541, A35
  28. Dedner, A., Kemm, F., Kröner, D., et al. 2002, Journal of Computational Physics, 175, 645
  29. Desch, S. J., & Mouschovias, T. C. 2001, ApJ, 550, 314
  30. Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269
  31. Duffin, D. F., & Pudritz, R. E. 2008, MNRAS, 391, 1659
  32. Ebert, R. 1955, Zeitschrift fur Astrophysik, 36, 222
  33. Enoch, M. L., Lee, J.-E., Harvey, P., Dunham, M. M., & Schnee, S. 2010, ApJ, 722, L33
  34. Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585
  35. Fiedler, R. A., & Mouschovias, T. C. 1993, ApJ, 415, 680
  36. Graves, D. T., Trebotich, D., Miller, G. H., & Colella, P. 2008, Journal of Computational Physics, 227, 4797
  37. Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
  38. Heitsch, F., & Hartmann, L. 2014, MNRAS, 443, 230
  39. Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25
  40. Hirano, N., & Liu, F.-c. 2014, ApJ, 789, 50
  41. Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128
  42. Joos, M., Hennebelle, P., Ciardi, A., & Fromang, S. 2013, A&A, 554, A17
  43. Krasnopolsky, R., Li, Z.-Y., & Shang, H. 2011, ApJ, 733, 54
  44. Krasnopolsky, R., Li, Z.-Y., Shang, H., & Zhao, B. 2012, ApJ, 757, 77
  45. Krumholz, M. R., Crutcher, R. M., & Hull, C. L. H. 2013, ApJ, 767, L11
  46. Kunz, M. W., & Mouschovias, T. C. 2009, ApJ, 693, 1895
  47. —. 2010, MNRAS, 408, 322
  48. Larson, R. B. 1969, MNRAS, 145, 271
  49. Lazarian, A., & Vishniac, E. T. 1999, ApJ, 517, 700
  50. Levermore, C. D. 1984, Journal of Quantitative Spectroscopy and Radiative Transfer, 31, 149
  51. Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321
  52. Li, H.-b., Fang, M., Henning, T., & Kainulainen, J. 2013, MNRAS, 436, 3707
  53. Li, H.-b., Goodman, A., Sridharan, T. K., et al. 2014a, ArXiv e-prints, arXiv:1404.2024
  54. Li, Z.-Y., Banerjee, R., Pudritz, R. E., et al. 2014b, ArXiv e-prints, arXiv:1401.2219
  55. Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180
  56. Li, Z.-Y., & McKee, C. F. 1996, ApJ, 464, 373
  57. Mac Low, M.-M., Norman, M. L., Konigl, A., & Wardle, M. 1995, ApJ, 442, 726
  58. Machida, M. N., & Hosokawa, T. 2013, MNRAS, 431, 1719
  59. Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2008a, ApJ, 676, 1088
  60. Machida, M. N., Inutsuka, S.-I., & Matsumoto, T. 2011a, PASJ, 63, 555
  61. Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2011b, ApJ, 729, 42
  62. —. 2014, MNRAS, 438, 2278
  63. Machida, M. N., & Matsumoto, T. 2011, MNRAS, 413, 2767
  64. Machida, M. N., Matsumoto, T., Hanawa, T., & Tomisaka, K. 2005a, MNRAS, 362, 382
  65. Machida, M. N., Matsumoto, T., Tomisaka, K., & Hanawa, T. 2005b, MNRAS, 362, 369
  66. Machida, M. N., Tomisaka, K., Matsumoto, T., & Inutsuka, S.-i. 2008b, ApJ, 677, 327
  67. Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24
  68. Masunaga, H., & Inutsuka, S.-i. 2000, ApJ, 531, 350
  69. Matsumoto, T., & Hanawa, T. 2003, ApJ, 583, 296
  70. Maury, A. J., André, P., Hennebelle, P., et al. 2010, A&A, 512, A40
  71. McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36
  72. Mellon, R. R., & Li, Z.-Y. 2008, ApJ, 681, 1356
  73. —. 2009, ApJ, 698, 922
  74. Mestel, L., & Spitzer, Jr., L. 1956, MNRAS, 116, 503
  75. Miyoshi, T., & Kusano, K. 2005, Journal of Computational Physics, 208, 315
  76. Mouschovias, T. C., & Paleologou, E. V. 1979, ApJ, 230, 204
  77. —. 1980, ApJ, 237, 877
  78. Mouschovias, T. C., & Spitzer, Jr., L. 1976, ApJ, 210, 326
  79. Nakano, T., Nishi, R., & Umebayashi, T. 2002, ApJ, 573, 199
  80. Ohashi, N., Saigo, K., Aso, Y., et al. 2014, ArXiv e-prints, arXiv:1410.0172
  81. Okuzumi, S. 2009, ApJ, 698, 1122
  82. Padovani, M., Galli, D., Hennebelle, P., Commerçon, B., & Joos, M. 2014, A&A, 571, A33
  83. Paleologou, E. V., & Mouschovias, T. C. 1983, ApJ, 275, 838
  84. Pezzuto, S., Elia, D., Schisano, E., et al. 2012, A&A, 547, A54
  85. Pineda, J. E., Arce, H. G., Schnee, S., et al. 2011, ApJ, 743, 201
  86. Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1
  87. Saigo, K., & Tomisaka, K. 2006, ApJ, 645, 381
  88. —. 2011, ApJ, 728, 78
  89. Saigo, K., Tomisaka, K., & Matsumoto, T. 2008, ApJ, 674, 997
  90. Sakai, N., Sakai, T., Hirota, T., et al. 2014, Nature, 507, 78
  91. Santos-Lima, R., de Gouveia Dal Pino, E. M., & Lazarian, A. 2012, ApJ, 747, 21
  92. —. 2013, MNRAS, 429, 3371
  93. Seaton, M. J., Yan, Y., Mihalas, D., & Pradhan, A. K. 1994, MNRAS, 266, 805
  94. Seifried, D., Banerjee, R., Pudritz, R. E., & Klessen, R. S. 2013, MNRAS, 432, 3320
  95. —. 2014, ArXiv e-prints, arXiv:1408.2989
  96. Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611
  97. Spruit, H. C., Stehle, R., & Papaloizou, J. C. B. 1995, MNRAS, 275, 1223
  98. Stamatellos, D., & Whitworth, A. P. 2009, MNRAS, 400, 1563
  99. Stamatellos, D., Whitworth, A. P., & Hubber, D. A. 2011, ApJ, 730, 32
  100. Stehle, R., & Spruit, H. C. 2001, MNRAS, 323, 587
  101. Susa, H., Doi, K., & Omukai, K. 2015, ArXiv e-prints, arXiv:1501.00087
  102. Tassis, K., & Mouschovias, T. C. 2005a, ApJ, 618, 769
  103. —. 2005b, ApJ, 618, 783
  104. —. 2007, ApJ, 660, 388
  105. Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2012, Nature, 492, 83
  106. —. 2013, ApJ, 771, 48
  107. Tokuda, K., Onishi, T., Saigo, K., et al. 2014, ApJ, 789, L4
  108. Tomida, K. 2014, ApJ, 786, 98
  109. Tomida, K., Machida, M. N., Saigo, K., Tomisaka, K., & Matsumoto, T. 2010a, ApJ, 725, L239
  110. Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2013, ApJ, 763, 6
  111. —. 2010b, ApJ, 714, L58
  112. Tomisaka, K. 1998, ApJ, 502, L163
  113. —. 2000, ApJ, 528, L41
  114. —. 2002, ApJ, 575, 306
  115. Tomisaka, K., & Tomida, K. 2011, PASJ, 63, 1151
  116. Toomre, A. 1964, ApJ, 139, 1217
  117. Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179
  118. Tsukamoto, Y., & Machida, M. N. 2011, MNRAS, 416, 591
  119. Umebayashi, T., & Nakano, T. 2009, ApJ, 690, 69
  120. Vorobyov, E. I., & Basu, S. 2006, ApJ, 650, 956
  121. —. 2010, ApJ, 719, 1896
  122. Wardle, M. 1991, MNRAS, 251, 119
  123. Wardle, M., & Koenigl, A. 1993, ApJ, 410, 218
  124. Wardle, M., & Ng, C. 1999, MNRAS, 303, 239
  125. Yorke, H. W., & Kaisig, M. 1995, Computer Physics Communications, 89, 29
  126. Ziegler, U., & Yorke, H. W. 1997, Computer Physics Communications, 101, 54
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description