RADIATION MAGNETOHYDRODYNAMIC SIMULATIONS OF PROTOSTELLAR COLLAPSE: PROTOSTELLAR CORE FORMATION

As stars are the most fundamental elements in the universe, their formation is one of the longstanding key topics in astrophysics. Because of its complicated nature and observational difficulties, computational simulations have played crucial roles in this field. Larson (1969) first showed the scenario from molecular cloud cores to protostellar cores with one-dimensional (1D) hydrodynamic simulations using the diffusion approximation for radiation transfer and found that two remarkable quasi-hydrostatic objects are formed in the course of protostellar collapse. A molecular cloud core initially collapses almost isothermally because thermal emission from dust grains is efficient in this phase. The gas temperature starts to rise when the central region gets dense and radiation cooling becomes inefficient, then the collapse almost stops and a quasi-hydrostatic object forms. This object, a so-called first (hydrostatic) core, evolves through accretion from the envelope. When the central temperature reaches about 2000 K where hydrogen molecules start to dissociate, the core becomes unstable and collapses again because H₂ dissociation is strongly endothermic and its gas pressure fails to balance gravity. This second collapse ends when molecular hydrogen dissociates completely and another quasi-static object is formed. This is the second core or the protostellar core. It acquires its mass via accretion on the dynamical timescale of the natal cloud core and evolves into a protostar and eventually a main-sequence star. More sophisticated simulations have been performed (Winkler & Newman 1980a, 1980b; Stahler et al. 1980a, 1980b, 1981; Masunaga & Inutsuka 2000) and the scenario of protostellar collapse under spherical symmetry is by now well established.

However, there is non-negligible internal motion (i.e., rotation and turbulence) in molecular cloud cores which makes the collapse completely different from the spherically symmetric models. Because centrifugal force is typically strong enough to prevent the formation of stars ("the angular momentum problem"), there must be sufficiently effective mechanisms of angular momentum transport. Gravitational torque via non-axisymmetric structure like spiral arms takes this role if magnetic fields are very weak. Bate (1998) first performed three-dimensional (3D) smoothed particle hydrodynamics (SPH) simulations and showed that gravitational torque redistributes angular momentum sufficiently and enables formation of protostellar cores. On the other hand, Tomisaka (1998, 2000, 2002) showed that magnetic fields transport angular momentum far more efficiently via torsional Alfvén waves (so-called magnetic braking; Mouschovias & Paleologou 1979, 1980) and bipolar outflows using two-dimensional (2D) MHD simulations. Machida et al. (2005a, 2005b) performed 3D nested-grid MHD simulations of similar situations and investigated the condition for fragmentation in the early phase before the second collapse begins. Machida et al. (2006a) investigated the evolution after the second collapse and showed that two different outflows are launched from the first core and the protostellar core. More detailed 3D simulations on this problem have been performed (Banerjee & Pudritz 2006; Hennebelle & Fromang 2008; Hennebelle & Teyssier 2008; Machida et al. 2008b) and effects of many physical ingredients such as misaligned magnetic fields (Matsumoto & Tomisaka 2004; Machida et al. 2006b; Ciardi & Hennebelle 2010; Joos et al. 2012), metallicities (Machida 2008), and turbulence (Seifried et al. 2012; Santos-Lima et al. 2012) have been investigated. Even weak magnetic fields can transport angular momentum very efficiently and modify the structure of the accretion flow significantly, and observations suggest that molecular clouds are considerably magnetized (Heiles & Crutcher 2005; Girart et al. 2006; Falgarone et al. 2008; Troland & Crutcher 2008; Crutcher et al. 2009). Therefore, magnetic fields are likely to play dominant roles in typical present-day environments.

Non-ideal MHD effects are also important in protostellar collapse. Typical magnetic flux in a molecular cloud core is far larger compared to that in a formed star, which means that there must be significant redistribution of magnetic flux during collapse. Moreover, Mellon & Li (2008, 2009) claimed that rotationally supported circumstellar disks are difficult to form in magnetized molecular cloud cores because the angular momentum transport due to magnetic fields is too efficient ("the magnetic braking catastrophe"; see also Li et al. 2011). Non-ideal MHD effects such as Ohmic dissipation and ambipolar diffusion are supposed to be responsible in this process because the ionization degree in star-forming clouds is very low. Ohmic dissipation redistributes sufficient magnetic flux outward from the centrally condensed high-density region and enables formation of circumstellar disks in the vicinity of protostars (Dapp & Basu 2010; Inutsuka et al. 2010; Machida et al. 2011; Machida & Matsumoto 2011). Gravitational torque becomes important in transporting angular momentum where magnetic fields are significantly weakened (Machida et al. 2010).

In those previous multi-dimensional simulations, they adopted the barotropic approximation for gas thermodynamics in which the thermal evolution is approximated with simple polytropic relations mimicking the results of 1D RHD simulations (e.g., Masunaga & Inutsuka 2000), instead of solving radiation transfer. Recently, Bate (2010, 2011; see also Whitehouse & Bate 2006) reported 3D SPH RHD simulations of formation of protostellar cores. Schönke & Tscharnuter (2011) also performed similar calculations of protostellar collapse using 2D axisymmetric RHD simulations. They commonly found that interesting phenomena happen when the protostellar cores are formed; bipolar outflows are launched from the disk-like first cores via radiation heating (not radiation force) from the protostellar cores. Considering the energy released in the second collapse and subsequent accretion onto the protostellar core, the irradiation is sufficient to heat up the gas in the first core and launch the outflows (Bate 2011). The barotropic approximation cannot reproduce such outflows driven by radiation heating. The structure of the first core which depends on the angular momentum distribution seems to be important in this phenomenon, and in the extreme case, the first core is almost blown up (Schönke & Tscharnuter 2011). Such violent phenomena may affect the story of protostellar collapse, circumstellar disk formation, and evolution of protostars.

Thus, both magnetic fields and radiation transfer are critically important in star formation processes. Commerçon et al. (2010) and Tomida et al. (2010b) independently performed 3D radiation magnetohydrodynamic (RMHD) simulations of protostellar collapse, but both stopped their calculations before the onset of the second collapse and they adopted the ideal MHD approximation. In this work, we report the first 3D RMHD simulations of protostellar collapse from molecular cloud cores to protostellar cores with and without Ohmic dissipation. To follow the evolution after the second collapse, we adopt a realistic equation of state (EOS) which takes chemical reactions into account. The goal of this study is to describe the realistic evolution in the early phase of protostellar collapse (i.e., until protostellar cores are formed) involving both radiation transfer and magnetic fields.

The plan of this paper is as follows. We describe the basic equations in Section 2. The adopted numerical methods are given in Section 3 and the simulation setups in Section 4. We show the results of our 3D RMHD simulations in Section 5. Section 6 is devoted to conclusions and discussions. In the Appendices, we explain the microphysical processes considered in our simulations.

In this work, we solve system equations including MHD, self-gravity, radiation transfer, a realistic EOS, and Ohmic dissipation. We use the comoving-frame RHD equations in which all the radiation related variables are defined in the comoving frame of the local fluid element (Castor 2004). For simplicity and to reduce computational load, we adopt the gray approximation (using frequency-averaged equations assuming the local thermodynamic equilibrium) and the flux-limited diffusion (FLD) approximation (Levermore & Pomraning 1981; Levermore 1984) for radiation transfer. We take Ohmic dissipation into account as the first step because it is supposed to be the most significant in the context of magnetic flux loss from high-density gas (Machida et al. 2007), although other non-ideal MHD effects such as ambipolar diffusion and Hall effects are also important (Nakano et al. 2002).

The governing equations are as follows:

$$\begin{eqnarray} &&\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf {v})=0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} && \frac{\partial \rho \mathbf {v}}{\partial t}+\nabla \cdot \left[\rho \mathbf {v}\otimes \mathbf {v}+\left(p+\frac{1}{2}|\mathbf {B}|^2\right)\mathbb {I}-\mathbf {B}\otimes \mathbf {B}\right]\nonumber\\ &&\quad =-\rho \nabla \Phi +\frac{\sigma _R}{c}\mathbf {F}_r, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&\frac{\partial \mathbf {B}}{\partial t}-\nabla \times (\mathbf {v}\times \mathbf {B}-\eta \nabla \times \mathbf {B})=0, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&\nabla \cdot \mathbf {B}=0, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&\frac{\partial e}{\partial t}+\nabla \cdot \left[\left(e+p+\frac{1}{2}|\mathbf {B}|^2\right)\mathbf {v}-\mathbf {B}(\mathbf {v}\cdot \mathbf {B})-\eta \mathbf {B}\times (\nabla \times \mathbf {B})\right]\nonumber \\ &&\quad =-\rho \mathbf {v}\cdot \nabla \Phi -c\sigma _P \big(a_rT_g^4-E_r \big)+\frac{\sigma _R}{c}\mathbf {F}_r\cdot \mathbf {v}, \end{eqnarray} \tag{ 5 }$$

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

$$\begin{eqnarray} &&\hspace{-30pt} \frac{\partial E_r}{\partial t}+\nabla \cdot [\mathbf {v}E_r]+\nabla \cdot \mathbf {F}_r+\mathbb {P}_r:\nabla \mathbf {v}=c\sigma _P\big(a_r T_g^4 -E_r\big), \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray*} &&\mathbf {F}_r=\frac{c\lambda }{\sigma _R}\nabla E_r,\hspace{10.0pt} \lambda (R)=\frac{2+R}{6+2R+R^2},\hspace{10.0pt} R=\frac{|\nabla E_r|}{\sigma _R E_r},\nonumber \\ &&\mathbb {P}_r=\mathbb {D}E_r,\hspace{10.0pt} \mathbb {D}=\frac{1-\chi }{2}\mathbb {I}+\frac{3\chi -1}{2}\mathbf {n}\otimes \mathbf {n},\hspace{10.0pt} \chi =\lambda +\lambda ^2R^2,\\ &&\quad \mathbf {n}=\frac{\nabla E_r}{|\nabla E_r|}\nonumber, \end{eqnarray*}$$

where ρ denotes the gas density, v the fluid velocity, B the magnetic flux density, p the gas pressure, T_g the gas temperature, $e=e_g+\frac{1}{2}\rho v^2+\frac{1}{2}|\mathbf {B}|^2$ the total gas energy density (e_g is the internal energy density of the gas), E_r the radiation energy density, Φ the gravitational potential, F_r the radiation energy flux, and $\mathbb {P}_r$ the radiation pressure tensor, respectively. c = 2.99792 × 10¹⁰ cm s⁻¹ is the speed of light and G = 6.673 × 10⁻⁸ cm³ g⁻¹ s⁻² is the gravitational constant. a_r = 4σ/c = 7.5657 × 10⁻¹⁵ erg cm⁻³ K⁻⁴ is the radiation (density) constant where σ = 5.6704 × 10⁻⁵ erg cm⁻² K⁻⁴ is the Stefan–Boltzmann constant. η is the resistivity and σ_R (σ_P) is the Rosseland (Planck) mean opacity. $\mathbb {I}$ denotes the unit matrix. ":" means the double dot product of two tensors, $\mathbb {A}:\mathbb {B}=\sum _i\sum _j A_{ij}B_{ji}$. These equations represent conservation of mass, the equation of motion, the induction equation including Ohmic dissipation, the solenoidal constraint, the Poisson's equation of gravity, and the gray FLD radiation transfer equations from top to bottom. Basically we use the Gaussian cgs units but we rescale the magnetic flux density to eliminate the constant coefficients, i.e., $\mathbf {B}=\mathbf {B}_0/\sqrt{4\pi }$ where B₀ is given in Gauss. Additionally, we need an EOS which gives the relations between the thermodynamic variables ρ, p, T, and e_g to close the system. We adopt the tabulated EOS in which we consider internal degrees of freedom and chemical reactions of seven species (H₂, H, H⁺, He, He⁺, He^{2 +}, and e⁻). Here, we assume the solar abundance, X = 0.7 and Y = 0.28. We also adopt the tabulated tables of the Rosseland and Planck mean opacities. In order to cover the huge dynamic range, we combine three opacity tables: Semenov et al. (2003), Seaton et al. (1994), and Ferguson et al. (2005). The resistivity is adopted from Umebayashi & Nakano (2009), while we additionally consider thermal ionization of potassium (K). The detailed descriptions about these microphysical processes are given in the Appendices.

Note that Machida et al. (2006a, 2007, 2008b) simplified the resistive term in Equation (3) to $\eta \nabla ^2\mathbf {B}$, but it is inadequate because the resistivity η has large spatial gradient according to the local gas density and temperature gradients, and what is worse, it violates the divergence free condition (Equation (4)). In this work, we correctly discretize the original equation and also include the effect of the resistivity in the energy equation (Equation (5)).

3.1. Operator Splitting

In order to solve the complex system described in the previous section, we divide the system into five parts and solve them separately on the nested-grid hierarchy. We update our system in the following strategy.

Step 1. Ideal MHD part:

$$\begin{eqnarray} \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf {v})&=&0, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \frac{\partial \rho \mathbf {v}}{\partial t}+\nabla \cdot \left[\rho \mathbf {v}\otimes \mathbf {v}+\left(p+\frac{1}{2}|\mathbf {B}|^2\right)\mathbb {I}-\mathbf {B}\otimes \mathbf {B}\right]&=&0, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \frac{\partial \mathbf {B}}{\partial t}-\nabla \times (\mathbf {v}\times \mathbf {B})&=&0, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \nabla \cdot \mathbf {B}&=&0, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \frac{\partial e}{\partial t}+\nabla \cdot \left[\left(e+p+\frac{1}{2}|\mathbf {B}|^2\right)\mathbf {v}-\mathbf {B}(\mathbf {v}\cdot \mathbf {B})\right]&=&0, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \frac{\partial E_r}{\partial t}+\nabla \cdot [\mathbf {v}E_r]&=&0. \end{eqnarray} \tag{ 13 }$$

Step 2. Self-gravity part:

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

$$\begin{eqnarray} \frac{\partial \rho \mathbf {v}}{\partial t}&=&-\rho \nabla \Phi, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} \frac{\partial e}{\partial t}&=&-\rho \mathbf {v}\cdot \nabla \Phi, \end{eqnarray} \tag{ 16 }$$

Step 3. Resistivity part:

$$\begin{eqnarray} \frac{\partial \mathbf {B}}{\partial t}+\nabla \times (\eta \nabla \times \mathbf {B})&=&0, \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} \frac{\partial e}{\partial t}-\nabla \cdot \left[\eta \mathbf {B}\times (\nabla \times \mathbf {B})\right]&=&0. \end{eqnarray} \tag{ 18 }$$

Step 4. Radiation part:

$$\begin{eqnarray} \frac{\partial e_g}{\partial t}&=&-c\sigma _P\big(a_rT_g^4-E_r\big), \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} \frac{\partial E_r}{\partial t}+\nabla \cdot \mathbf {F}_r+\mathbb {P}_r:\nabla \mathbf {v}&=&c\sigma _P \big(a_r T_g^4 - E_r\big), \end{eqnarray} \tag{ 20 }$$

Step 5. Radiation force part:

$$\begin{eqnarray} \frac{\partial \rho \mathbf {v}}{\partial t}&=&\frac{\sigma _R}{c}\mathbf {F}_r, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \frac{\partial e}{\partial t}&=&\frac{\sigma _R}{c}\mathbf {F}_r\cdot \mathbf {v}. \end{eqnarray} \tag{ 22 }$$

Introducing such an operator-splitting technique causes loss of time accuracy, and the overall accuracy of our time-integration scheme is first order.

3.2. Magnetohydrodynamics and Self-Gravity

3.3. Ohmic Dissipation

Because of the low ionization degree, non-ideal MHD effects such as Ohmic dissipation, ambipolar diffusion, and the Hall effect work during protostellar collapse. Ohmic dissipation becomes dominant in the high-density region (typically ρ ≳ 10⁻¹¹ g cm⁻³) and ambipolar diffusion dominates in the low-density region (Nakano et al. 2002; Kunz & Mouschovias 2009, 2010), but in this work we focus on Ohmic dissipation as the first step because it is supposed to be most significant in the context of the magnetic flux problem and the angular momentum problem (Machida et al. 2011; Dapp et al. 2012). The introduction of ambipolar diffusion will be a future work. The timescale of Ohmic dissipation becomes shorter than the dynamical timescale in the high-density region (ρ ≳ 10⁻¹⁰ g cm⁻³; Figure 28) and significant diffusion of magnetic flux occurs (Nakano et al. 2002). Note that these non-ideal MHD processes do not reduce total magnetic flux but redistribute centrally condensed magnetic flux outward from the high-density region (see below, Figure 16), although this process is often mentioned as "magnetic flux loss."

We discretize (17) and (18) in the same way as Matsumoto (2011). Additionally, we introduce the correction terms proportional to ∇ · B to improve the nature of the system as described in Graves et al. (2008). Since the timescale of the diffusion of magnetic fields can be far shorter than the dynamical timescale (Figure 28), here we adopt the Super-Time-Stepping (STS) method proposed by Alexiades et al. (1996; for astrophysical applications, see O'Sullivan & Downes 2006; Choi et al. 2009; Commerçon et al. 2011b) when the timestep for the resistivity part is shorter than the hydrodynamic timestep. STS is constructed based on the simple explicit scheme, but we can relax the Courant–Friedrich–Levy (CFL) stability condition by using the following sequence of sub-timesteps:

$$\begin{eqnarray} &&\tau _j=\Delta t_{{\rm exp}}\left[(-1+\nu)\cos \left(\frac{2j-1}{N}\frac{\pi }{2}\right)+1+\nu \right]^{-1}. \end{eqnarray} \tag{ 23 }$$

The longer time integration of ΔT_STS is achieved after N sub-timesteps,

$$\begin{eqnarray} &&\Delta T_{{\rm STS}}=\sum _j^{N}\tau _j=\Delta t_{{\rm exp}}\frac{N}{2\sqrt{\nu }}\left[\frac{(1+\sqrt{\nu })^{2N}-(1-\sqrt{\nu })^{2N}}{(1+\sqrt{\nu })^{2N}+(1-\sqrt{\nu })^{2N}}\right],\nonumber\\ \end{eqnarray} \tag{ 24 }$$

where Δt_exp is the explicit timestep that satisfies the CFL condition and ν is a small positive parameter that controls the stability and efficiency of the scheme. Smaller ν gives better acceleration but STS may give unphysical results when it is too small. The optimal choice of ν depends on the problem, but typically ν ∼ 0.01 seems to be good for parabolic problems. From Equation (24), the maximum acceleration for given ν compared to the explicit scheme (ΔT_exp = NΔt_exp) is $0.5/\sqrt{\nu }$ and the most efficient calculation is realized around $N\sim 0.5/\sqrt{\nu }$. In this work, we adopt ν = 0.01 and N = 6. Even in the most time-consuming situation (i.e., the magnetic Reynolds number is very low, ρ ∼ 10⁻⁹ g cm⁻³) in our production runs, STS keeps the computational cost lower than or similar to the total of all the other parts (including MHD, self-gravity, and radiation transfer) at most.

3.4. Radiation Transfer

The timescale related to radiation can be far shorter than that related to hydrodynamics in star formation simulations. The light speed c = 2.99792 × 10¹⁰ cm s⁻¹ is far larger than typical fluid velocity in star formation, which is on the order of (1 – 100) × 10⁵ cm s⁻¹. Moreover, the radiation source terms are strongly nonlinear. Therefore, the differential equations of the radiation subsystem can be very stiff. It is unreasonable to calculate such a system using explicit schemes or even STS. Therefore, we adopt an implicit time-integration scheme that is stable regardless of the timescale of involved physical processes.

3.4.1. Discretization

We adopt the first-order backward Euler method that is simple and stable. Here, we discretize (19) and (20) in one dimension (extension to multi-dimension is straightforward) as follows:

$$\begin{eqnarray} \frac{e_{g,i}^{n+1}-e_{g,i}^n}{\Delta t}&=&-c\sigma _P^n \big[a_r \big(T_{g,i}^{n+1}\big)^4-E_{r,i}^{n+1}\big], \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} \frac{E_{r,i}^{n+1}-E_{r,i}^n}{\Delta t}&-&\frac{1}{\Delta x}\bigg[\bigg(\frac{c\lambda }{\sigma _R}\bigg)_{i+\frac{1}{2}}^n \frac{E_{r,i+1}^{n+1}-E_{r,i}^{n+1}}{\Delta x}\nonumber\\ \hphantom{\frac{E_{r,i}^{n+1}-E_{r,i}^n}{\Delta t}-}&-&\bigg(\frac{c\lambda }{\sigma _R}\bigg)_{i-\frac{1}{2}}^n \frac{E_{r,i}^{n+1}-E_{r,i-1}^{n+1}}{\Delta x}\bigg]\nonumber \\ \hphantom{\frac{E_{r,i}^{n+1}-E_{r,i}^n}{\Delta t}}&=&c\sigma _P^n \big[a_r \big(T_{g,i}^{n+1}\big)^4-E_{r,i}^{n+1}\big]-\mathbb {P}_{r,i}^{n+1}:(\nabla \mathbf {v})_i.\qquad \end{eqnarray} \tag{ 26 }$$

Superscripts and subscripts denote the indices of the discretized time and space, respectively. To avoid numerical difficulties due to strong nonlinearity in the flux limiter and the opacities and to achieve stable integration, we adopt the time-lagged opacities and flux limiter (Castor 2004). This may cause a loss of the time accuracy of the scheme, but we confirmed that it does not matter in our production runs because the timescale of the evolution of radiation fields is similar to the hydrodynamic timescale and therefore well resolved. We also use the time-lagged Eddington Tensor $\mathbb {D}^n$, which yields the radiation energy tensor $\mathbb {P}_r^{n+1}=\mathbb {D}^n E_r^{n+1}$.

There are some ways to evaluate the opacities and the flux limiter at the cell interfaces, i + 1/2. Here, we follow Howell & Greenough (2003) and adopt the surface formula that gives good flux even at a sharp surface of optically thick material like the surface of a first core:

$$\begin{eqnarray} &&\hspace{-4pt}\sigma _{R,i+\frac{1}{2}}=\min \left[\frac{\sigma _{R,i}+\sigma _{R,i+1}}{2},\max {\left(\frac{2\sigma _{R,i}\sigma _{R,i+1}}{\sigma _{R,i}+\sigma _{R,i+1}},\frac{4}{3\Delta x}\right)}\right].\nonumber\\ \end{eqnarray} \tag{ 27 }$$

We evaluate the flux limiter at the cell interface $\lambda (R)_{i+\frac{1}{2}}$ using the following equations:

$$\begin{eqnarray} &&R_{i+\frac{1}{2}}=\frac{|(\nabla E_r)_{i+\frac{1}{2}}|}{\sigma _{R,i+\frac{1}{2}}{E_{r,i+\frac{1}{2}}}}, \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} &&E_{r,i+\frac{1}{2}}=\frac{E_{r,i}+E_{r,i+1}}{2}, \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&(\nabla E_r)_{i+\frac{1}{2}}=\frac{E_{r,i+1}-E_{r,i}}{\Delta x}. \end{eqnarray} \tag{ 30 }$$

In three dimensions, the radiation energy gradient should be replaced with

$$\begin{eqnarray} \hspace{-16pt}&& (\nabla E_r)_{i+\frac{1}{2},j,k}\nonumber\\ \hspace{-14pt}&&\quad =\left(\begin{array}{@{}c@{}}\displaystyle \frac{E_{r,i+1}-E_{r,i}}{\Delta x}\\ \displaystyle \frac{E_{r,i+1,j+1,k}+E_{r,i,j+1,k}-E_{r,i+1,j-1,k}-E_{r,i,j-1,k}}{4\Delta y}\\ \displaystyle \frac{E_{r,i+1,j,k+1}+E_{r,i,j,k+1}-E_{r,i+1,j,k-1}-E_{r,i,j,k-1}}{4\Delta z} \end{array} \right)\!.\nonumber\\ \end{eqnarray} \tag{ 31 }$$

3.4.2. Newton–Raphson Iterations

To solve the nonlinear system (25) and (26), we perform the Newton–Raphson iterations (Press et al. 2007). In this method, we search for the zero-points of the residual functions f_i(X). We can find the root iteratively using the following matrix equation based on the Taylor expansion:

$$\begin{eqnarray} &&\sum _{j=1}^N\frac{\partial f_i}{\partial x_j}\delta X_j=-f_i(\mathbf {X}). \end{eqnarray} \tag{ 32 }$$

In our system, X is the vector of the gas and radiation energy densities in all the cells:

$$\begin{eqnarray} &&\hspace{-16pt}\mathbf {X}=(e_{g,1},\ldots,e_{g,i},e_{g,i+1},\ldots,E_{r,1},\ldots,E_{r,i},E_{r,i+1},\ldots)^T.\nonumber\\ \end{eqnarray} \tag{ 33 }$$

The residual functions are given from (25) and (26) as follows:

$$\begin{eqnarray} &&f^g_i=e_{g,i}^{n+1}-e_{g,i}^n+\Delta t c\sigma _P^n \big[a_r \big\lbrace T_g \big(\rho _i,e_{g,i}^{n+1} \big)\big\rbrace ^4-E_{r,i}^{n+1}\big], \quad \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} f^r_i&=&E_{r,i}^{n+1}-E_{r,i}^n+\Delta t\Biggl [\frac{1}{\Delta x}\Biggl\lbrace \Biggl(\frac{c\lambda }{\sigma _R}\Biggr)_{i+\frac{1}{2}}^n\frac{E_{r,i+1}^{n+1}-E_{r,i}^{n+1}}{\Delta x}\nonumber\\ &&-\Biggl(\frac{c\lambda }{\sigma _R}\Biggr)_{i-\frac{1}{2}}^n\frac{E_{r,i}^{n+1}-E_{r,i-1}^{n+1}}{\Delta x}\Biggr\rbrace \nonumber \\ &&-c\sigma _P^n \big[a_r \big\lbrace T_g \big(\rho _i,e_{g,i}^{n+1} \big)\big\rbrace ^4-E_{r,i}^{n+1}\big]+\{ \mathbb {D}^n_i:(\nabla \mathbf {v})_i\}E_{r,i}^{n+1}\big ].\nonumber\\ \end{eqnarray} \tag{ 35 }$$

We can rewrite (32) explicitly:

$$\begin{eqnarray} &&\frac{\partial f_i^g}{\partial e_{g,i}^{n+1}}\delta e_{g,i}^{n+1}+\frac{\partial f_i^g}{\partial E_{r,i}^{n+1}}\delta E_{r,i}^{n+1}=-f_i^g, \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} && \frac{\partial f_i^r}{\partial e_{g,i}^{n+1}}\delta e_{g,i}^{n+1}+\frac{\partial f_i^r}{\partial E_{r,i}^{n+1}}\delta E_{r,i}^{n+1}+\frac{\partial f_i^r}{\partial E_{r,i+1}^{n+1}}\delta E_{r,i+1}^{n+1}\nonumber\\ &&\quad +\frac{\partial f_i^r}{\partial E_{r,i-1}^{n+1}}\delta E_{r,i-1}^{n+1}=-f_i^r. \end{eqnarray} \tag{ 37 }$$

By substituting (36) into (37), we can eliminate the equation related to gas energy (Hayes et al. 2006; this procedure corresponds to performing partial LU decomposition analytically). Then we obtain the matrix equation:

$$\begin{eqnarray} &&\left(\frac{\partial f_i^r}{\partial E_{r,i}^{n+1}}-\frac{\partial f_i^r}{\partial e_{g,i}^{n+1}}\frac{\partial f_i^g}{\partial E_{r,i}^{n+1}}\bigg {/}\frac{\partial f_i^g}{\partial e_{g,i}^{n+1}}\right)\delta E_{r,i}^{n+1}+\frac{\partial f_i^r}{\partial E_{r,i+1}^{n+1}}\delta E_{r,i+1}^{n+1}\nonumber\\ &&\quad+ \frac{\partial f_i^r}{\partial E_{r,i-1}^{n+1}}\delta E_{r,i-1}^{n+1} =\frac{\partial f_i^r}{\partial e_{g,i}^{n+1}}\!\Bigg {/}\frac{\partial f_i^g}{\partial e_{g,i}^{n+1}}f_i^g-f_i^r. \end{eqnarray} \tag{ 38 }$$

The derivatives are given as follows:

$$\begin{eqnarray} \hspace{-20pt}\frac{\partial f_i^g}{\partial e_{g,i}^{n+1}}&=&1+4\Delta t\, c\sigma _P^n a_r \big[T_g\big(\rho _i,e_{g,i}^{n+1}\big)\big]^3\frac{\partial T_g}{\partial e_{g}}\big(\rho _i,e_{g,i}^{n+1}\big), \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} \frac{\partial f_i^g}{\partial E_{r,i}^{n+1}}&=&-\Delta t\, c \sigma^n _P, \end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray} \hspace{-16pt}\frac{\partial f_i^r}{\partial e_{g,i}^{n+1}}&=&{-}4\Delta t\, c\sigma^n _P a_r \big[T_g\big(\rho _i,e_{g,i}^{n+1}\big)\big]^3\frac{\partial T_g}{\partial e_{g}} \big(\rho _i,e_{g,i}^{n+1}\big), \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} \hspace{-14pt}\frac{\partial f_i^r}{\partial E_{r,i}^{n+1}}&=&1+\Delta t \Biggl[c \sigma^n _P +\mathbb {D}_i^n:(\nabla \mathbf {v})_i \nonumber\\ &&-\frac{c}{\Delta x^2}\left\lbrace \left(\frac{\lambda }{\sigma _R}\right)_{i+\frac{1}{2}}^n+\left(\frac{\lambda }{\sigma _R}\right)_{i-\frac{1}{2}}^n\right\rbrace \Biggr], \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} \frac{\partial f_i^r}{\partial E_{r,i\pm 1}^{n+1}}&=&\frac{c\Delta t}{\Delta x^2}\left(\frac{\lambda }{\sigma _R}\right)_{i\pm \frac{1}{2}}^n. \end{eqnarray} \tag{ 43 }$$

T_g(ρ, e_g) and (∂T_g/∂e_g)(ρ, e_g) are given from the tabulated EOS. Note that this Jacobi matrix is symmetric. We can obtain the solution by updating E_{r, i} and e_{g, i} using δE_{r, i} and δe_{g, i} calculated from Equations (38) and (36) until f_i and δX become sufficiently small. As the initial guess for E^{n + 1}_{r, i} and e^{n + 1}_{g, i}, we adopt the solution at the timestep n. In our simulations, we use convergence thresholds like max (f^g_i/e_{g, i}, f^r_i/E_{r, i}) < 5 × 10⁻⁴ and max (δe_{g, i}/e_{g, i}, δE_{r, i}/E_{r, i}) < 5 × 10⁻⁴. If these thresholds are not satisfied after many iterations, we take substeps with shorter Δt and try again until we obtain the converged solution successfully.

3.4.3. Linear System Solver

In 3D Cartesian coordinate, the Jacobi matrix in (38) is a very large sparse seven-diagonal matrix. Moreover, because of the strong nonlinearity of the system, the matrix is not necessarily diagonally dominant. Therefore, we need an efficient and robust sparse matrix solver and extensive computational resources to solve this large (typically 64³ = 262144 cells per grid level) linear system. We found that the combination of the BiCGStab (stabilized bi-conjugate gradient) solver and the incomplete LU (ILU) decomposition preconditioner without fill-in works well for our problems.⁷

3.4.4. Radiation Force

We simply integrate the radiation force terms in Step 5 using the obtained solution in Step 4. These terms are relatively small, at least in the early phase of low-mass star formation processes.

3.5. Nested-Grid

3.6. Boundary Conditions

We use unstable Bonnor–Ebert (hereafter BE; Bonnor 1956; Ebert 1955) spheres of T = 10 K as the initial conditions of our simulations, mimicking isolated molecular cloud cores. We construct a critical BE sphere (with a dimension less radius of ξ = 6.45) and make it unstable by increasing the gas density by a factor of A₀. Then we introduce uniform (rigid-body) rotation, uniform magnetic fields, and m = 2 density perturbation, where m is the number of longitudinal modes. The initial density profile is given as follows:

$$\begin{eqnarray} &&\rho (r)=\rho _c\rho _0(r)(1+A_0)\left[1+A_2\frac{r^2}{R^2}\cos (2\phi)\right], \end{eqnarray} \tag{ 44 }$$

where ρ_c is the gas density at the center of the cloud core, ρ₀(r) is the normalized density profile of the critical BE sphere, R is the radius of the critical BE sphere, and A₂ is the amplitude of m = 2 perturbation, respectively. In order to minimize the effect of initial resolution, we adopt this "regularized" m = 2 perturbation which is smooth at the center of the cloud in contrast to Boss & Bodenheimer (1979).

In this work, we adopt ρ_c = 10⁻¹⁸ g cm⁻³ and A₀ = 0.2, which yield an unstable BE-like sphere whose mass and radius are M ∼ 1 M_☉ and R ∼ 4.25 × 10⁻² pc ∼ 8800 AU. The initial free-fall time at the center of the cloud is t_ff ∼ 6.08 × 10⁴ yr.

We calculated a spherically symmetric model and four magnetized rotating models with and without Ohmic dissipation; the parameters of the models are summarized in Table 1. The first letter of the model name denotes the treatment of magnetic fields: I denotes an ideal MHD model and R a resistive MHD model. The second letter represents the initial rotation speed: F is fast and S is slow. Note that we choose these parameters so that the first core disks do not fragment before the second collapse because of the limitation of our nested-grids. Model SP is the spherical model without rotation and magnetic fields. For magnetized models, we impose the uniform magnetic fields of 20 μG parallel to the rotation axis. The corresponding mass-to-flux ratio normalized by the critical value of stability is μ₀ ≡ (M/Φ)/(M/Φ)_crit ∼ 3.8 where Φ = πR²B₀ and (M/Φ)_crit = (0.53/3π)(5/G)^1/2. Here we adopt the critical mass-to-flux ratio of Mouschovias & Spitzer (1976), but we should regard this value as just a guide because our initial conditions are not uniform. There is another similar threshold for stability between gravity and magnetization derived for disks (Nakano & Nakamura 1978); we have λ₀ ≡ (Σ/B)/(Σ/B)_crit = 7.2 at the center of the cloud, where (Σ/B)_crit = (4π²G)^−1/2. These mass-to-flux ratios indicate that our magnetized models are in the magnetically super-critical regime but considerably magnetized.

We define the origin of time as the epoch of formation of the protostellar core, i.e., when the central gas density exceeds ρ_c = 10⁻³ g cm⁻³ for the first time, for descriptive purposes. We stop our simulations when the central temperature reaches T_c ∼ 10⁵ K. Our simulations show the formation and earliest evolution of protostars. The typical spatial resolution around the surface of the first cores (r ∼ 3 − 5 AU) is Δx ∼ 0.14 AU (L = 12) and the finest resolution at the end of the simulations is Δx ∼ 6.6 × 10⁻⁵ AU ∼ 0.014 R_☉ (L = 23). The nested-grid technique enables us to achieve such a huge spatial dynamic range of more than eight orders of magnitude.

5.1. Spherical Model

We first show the results of the spherical model SP to understand the whole evolution from a molecular cloud core to a protostellar core, and to demonstrate validity of our code.

5.1.1. Thermal Evolution and EOS

Before discussing the global structure and evolution, we explain the thermal evolution of the central gas element. To see the effects of the microphysics involved in the EOS on the thermal evolution, we show the evolution track of the central gas element in the ρ − T plane and the adiabatic index Γ in Figure 1. The evolution is consistent with the results of the more elaborate 1D spherically symmetric RHD simulations (Masunaga & Inutsuka 2000) except for the details of the EOS. From this figure, we can understand the relation between the thermal evolution of the gas and the microphysical processes.

Zoom In Zoom Out Reset image size

While the gas density is low (ρ_c ≲ 10⁻¹³ g cm⁻³), the gas collapses isothermally because radiation cooling is very efficient, and the EOS does not matter in this regime. When the gas density gets sufficiently high and radiation cooling becomes inefficient, the gas evolves quasi-adiabatically and the temperature rises. Beyond this point, the EOS plays almost dominant roles in the thermal evolution of the gas (at the center of the cloud; note that the gas in outer region is affected by radiation transfer). While the gas is still cold (T_c ≲ 100 K), the adiabatic index Γ is about 5/3 because the rotational transitions cannot be excited and molecular hydrogen behaves like monoatomic molecules. Γ decreases to ∼7/5 when the gas becomes warm enough to excite rotation. When the temperature exceeds T_c ∼ 2000 K, molecular hydrogen starts to dissociate and the second collapse begins. The endothermic reaction of hydrogen molecule dissociation significantly decreases the effective adiabatic index to Γ ∼ 1.1, below the critical value of stability of a self-gravitational sphere, Γ_crit = 4/3. After the completion of the dissociation, the gas evolves quasi-adiabatically again, while the ionizations of hydrogen and helium slightly affect the evolution, making the adiabatic index softer. From Figure 1, we can see that the ionizations of hydrogen and helium do not cause the "third collapse"; those endothermic reactions proceed gradually because the reverse reactions (recombinations) occur rapidly in the high-density regions and the adiabatic index remains larger than the critical value.

In this work we assume the ortho:para ratio of molecular hydrogen to be 3:1, but if we adopt the equilibrium ratio, additional energy is consumed to convert parahydrogen to orthohydrogen (ΔE/k ∼ 170 K), resulting in the softer adiabatic index (Boley et al. 2007; Stamatellos & Whitworth 2009). The gas does not forget the thermal history in this phase because its evolution is almost adiabatic, therefore the ortho–para ratio has non-negligible impact on the whole star formation process. We must keep in mind that the differences in the EOS quantitatively affect the evolution and properties of the first and protostellar cores such as their radii, masses, and stability against fragmentation.

5.1.2. First Core

The first core is formed at t ∼ −650 yr when the central density exceeds ρ_c ∼ 10⁻¹² g cm⁻³. Its radius is initially about ∼5.4 AU, and it contracts gradually to ∼3 AU as it evolves (Figure 2). The evolution and structure of the first core are in good agreement with the results of the 1D spherically symmetric RHD simulations using full non-gray radiation transfer (Masunaga et al. 1998; Masunaga & Inutsuka 2000) and the 3D SPH RHD simulations using the gray FLD approximation (Whitehouse & Bate 2006). The outer region attains a higher entropy than the central gas element because of heating at the shock and radiation from the central hot region (Figure 3). The shock at the surface of the first core is isothermal, i.e., it is a super-critical shock at which almost all the kinetic energy is radiated away to the upstream, as discussed in Commerçon et al. (2011a). The lifetime of the spherical first core is a bit longer than that in Bate (2011), probably because of the different initial conditions; Bate (2011) adopted a uniform sphere with higher central gas density, which is more unstable and gives higher accretion rate than our BE sphere.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When the gas temperature exceeds the evaporation temperature of all the dust components (T ∼ 1500 K), the opacities around the central region drop significantly. The temperature distribution within this dust free region becomes almost flat. As Figure 2 indicates, the dust evaporation front is located at R ∼ 1.2 AU at the end of the first core phase (t ∼ −2.7 yr). The dust evaporation front is also visible as a small jump in Figure 3. Schönke & Tscharnuter (2011) pointed out the importance of the dust evaporation on the dynamics of first cores, but it is not prominent in the spherically symmetric model.

5.1.3. Protostellar Core

The second collapse begins when the central temperature exceeds T_c ∼ 2000 K. Soon after the onset of the second collapse, the protostellar core is formed in the short dynamical timescale of several years (the free-fall time corresponding to the central density when the second collapse starts (ρ_c ∼ 10⁻⁸ g cm⁻³) is only about 0.67 years). Within 0.7 years after the formation of the protostellar core, it acquires M_PC ∼ 2 × 10⁻² M_☉ and the averaged accretion rate is very high, 2.7 × 10⁻²M_☉ yr⁻¹. The protostellar core expands due to the addition of newly accreted gas and the accretion shock at the surface of the protostellar core quickly propagates outward. The jump condition at this shock is almost adiabatic, which means that the flow is radiatively inefficient ("hot accretion") in this early phase. The outer region of the protostellar core is heated up by the shock and attain a high entropy, leading the protostellar core to be convectively stable (Stahler et al. 1980a, 1980b). These thermal properties and expansion of protostellar cores cannot be reproduced by the barotropic approximation in which the shock heating is not taken into account. At the end of the simulation, the radius of the protostellar core is about R_PC ∼ 0.047 AU ∼ 10 R_☉. From the virial theorem, the energy released in this phase can be estimated to be (GM²_PC/2R_PC) ∼ 6.8 × 10⁴³ erg, which is consistent with the total dissociation energy of molecular hydrogen, (XM_PC/2m_H)χ_dis ∼ 5.7 × 10⁴³ erg where χ_dis = 7.17 × 10⁻¹² erg is the dissociation energy of H₂ (Liu et al. 2009).

This radius at the end of the simulation is about 2.5 times larger than the radius of the protostellar core obtained in Masunaga & Inutsuka (2000; ∼4 R_☉), and still increasing. The expansion in the adiabatic accretion phase had been already reported in Larson (1969) and also discussed in Stahler et al. (1986; see also Winkler & Newman 1980b), and our radius is consistent with their results. It continues to expand until almost all the gas in the first core has accreted onto the protostellar core when the accretion rate gets significantly low and the optical depth of the envelope becomes low enough for radiation cooling. It will take about the average free-fall time of the whole first core, t_{ff, FC} = (3π/32Gρ_FC)^1/2 ∼ 5.5 yr where ρ_FC = (3M_FC/4πR³_FC) = 1.5 × 10⁻¹⁰ g cm⁻³, M_FC ∼ 3 × 10⁻²M_☉, and R_FC ∼ 3 AU. Therefore, this expansion is a highly transient phenomenon. In other words, our protostellar core has not yet settled. Note that the gas does not flow outward in this expansion phase, but the newly accreted gas is loaded on top of the protostellar core. It is different from the violent explosion, or "hiccup" of protostellar cores discussed years ago (e.g., Boss 1989).

5.2. Rotating Models

5.2.1. Overview

We show the evolution tracks of the central gas density and temperature as functions of time in Figure 4, and the evolution tracks in the ρ − T plane in Figure 5. The first cores are formed when the central density exceeds ρ_c ≳ 10⁻¹² g cm⁻³ and the second collapse begins when ρ_c ≳ 10⁻⁸ g cm⁻³. The lifetimes of the first cores are about 650, 720, 800, 850, and 950 years in SP, IS, IF, RS, and RF, respectively. The presence of rotation extends the first core lifetime but its effect is not significant compared to non-magnetized cases (Saigo et al. 2008; Bate 2011; Tomida et al. 2010a). The resistive models have slightly longer lifetimes because magnetic fields are weakened by Ohmic dissipation and the efficiency of the angular momentum transport is reduced, but they are still not very long-lived.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Interestingly, all the models show essentially the same thermal evolution (Figure 5). This is because the central regions of the first cores contract almost spherically because of the efficient angular momentum transport via magnetic fields and their evolution converges on the Larson–Penston self-similar solution (Larson 1969; Penston 1969). In the central regions of the collapsing cloud cores, the thermal energy is dominantly larger than the magnetic energy and dissipation of magnetic fields (Joule heating) does not affect the thermal evolution significantly. Hence the thermal properties of the formed protostellar cores, such as gas entropy, have already been set in the first core phase to some extent and seem to depend weakly on some initial cloud parameters, such as rotation and magnetic fields, when efficient mechanisms of angular momentum transport are present. Without magnetic fields, the initial conditions have larger impact on the thermal evolutions, particularly when fragmentation occurs due to fast initial rotation (Bate 2011). The thermal evolution may also depend on many other factors such as initial temperature, dust opacities, and mass of the natal molecular cloud core (Tomida et al. 2010a).

5.2.2. Outflows and First Cores

After the formation of the first cores, slow and loosely collimated outflows are launched from the first cores by magneto-centrifugal force (Blandford & Payne 1982; Kudoh et al. 1998). We show the density and temperature distributions of the outflow scale (corresponding to the grid level l = 8 or ∼140 AU) and of the first core scale (l = 11 or ∼18 AU) in Figures 6–13 at the end of the first core phase. The outflows can be seen as slowly outgoing gas denser than the envelope that extends to z = 55, 70, 60, and 80 AU in IS, IF, RS, and RF, respectively. The outflows are not prominent in the temperature plots except for the thin shock-heated layers because they are optically thin and radiation cooling/heating is efficient. We also show the profiles of the gas density, temperature, and velocities along the x- and z-axes in Figure 14. Roughly speaking, the properties of the first cores and outflows are consistent with those in previous studies (Commerçon et al. 2010; Tomida et al. 2010b, etc.).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The properties of the outflows such as velocities and traveling distances are similar in all the models. This is because the driving radii are similar in all the models, ∼10 AU, where the gas density is not high enough for the resistivity to work. The outflow velocities are comparable to the rotational velocities at this radius, ∼1 km s⁻¹. Therefore, the traveling distances are almost proportional to the lifetimes of the first cores. It seems difficult to find the effect of the resistivity on the outflow properties because there is no essential difference in the outflows between the resistive and ideal MHD models (see also Yamada et al. 2009).

Although the outflows are similar, the first cores look quite different between the resistive and ideal MHD models. In the ideal models, the first cores are virtually non-rotating because of the efficient magnetic braking (R ≲ 3 AU in Figure 14). However, the magnetic braking is less efficient in the resistive models and there remains a considerable amount of angular momentum. To show the effects of the resistivity clearly, we plot the total angular momentum within the first cores in Figure 15. Here, we simply measure the angular momentum where the gas density is above a critical value, ρ_crit = 10⁻¹³ g cm⁻³. The resistivity is efficient where ρ ≳ 10⁻¹⁰ g cm⁻³. The magnetic flux is extracted from the central high-density region and redistributed to the outer thinner region, but the redistributed magnetic fields do not show any significant effect dynamically, at least in this phase. In both models with fast and slow rotation, the first cores in the resistive models attain about twice larger angular momenta compared with the corresponding ideal MHD models. Even though these differences do not have a significant impact on the evolution of the first cores, they become important later in the protostellar core phase.

Zoom In Zoom Out Reset image size

The outstanding difference between the ideal and resistive models is that the first cores and the surrounding (pseudo) disks warp in the ideal MHD models. These are likely to be the magnetic interchange instability (Spruit & Taam 1990; Spruit et al. 1995; Stehle & Spruit 2001; Krasnopolsky et al. 2012). Uniformly rotating disks are unstable when the ratio between the poloidal magnetic flux density and the surface density decreases radially; (d(B_z/Σ)/dr) < 0. Both ideal and resistive models satisfy this condition in 1 AU ≲ R ≲ 10 AU (Figure 16), though only the ideal MHD models suffer from the warping in our simulations. The timescale of growth is probably long in the resistive cases or we need to consider more realistic situations including the effects of finite disk thickness and gas infall to understand these phenomena, but they are beyond the scope of this paper. Possibly, the magnetically driven warping instability (Lai 2003) may also contribute to the warping. We should note that even small vertical perturbation induces artificial reconnection of magnetic fields at the mid-plane because the magnetic fields are pinched onto the disk mid-plane and the directions of the magnetic fields above and below the mid-plane are anti-parallel, then things become chaotic. Because this warping does not affect the evolution of the protostellar cores, we do not discuss it further in this work.

Zoom In Zoom Out Reset image size

We can see that the gas within the dust evaporation front slowly infalls due to loss of the pressure gradient in all the models. However, the dust evaporation front is still confined in the first core and therefore its effects do not seem significant. This is different from the non-magnetized RHD simulations where the front expands beyond the first core surface (Schönke & Tscharnuter 2011). In our simulations, the angular momentum transport is very efficient and the first core properties are quite similar to the spherically symmetric cases, even in the resistive models.

Interestingly, the size (height) of the first core is slightly larger in the resistive models. It is about 3 AU in the ideal MHD models, which is similar to the spherical model, but is about 5 AU in the resistive models. This is interpreted as a consequence of energy transport and additional heating by Ohmic dissipation. The magnetic fields are transported outward, then heat up and inflate the outer region (note that the resistivity is most effective around ρ ∼ 10⁻⁹ g cm⁻³).

5.2.3. Protostellar Cores and Jets

We show the density and temperature cross sections of the protostellar core scale (l = 15 or 16, corresponding to ∼1.1 AU or ∼0.54 AU) at the end of the simulations in Figures 17–20. We also show the profiles along the x- and z-axes in Figure 21. We stop our simulations when the central temperature reaches T_c ∼ 10⁵ K, corresponding to 1.05, 0.44, 0.90, and 1.25 years after the formation of the protostellar cores in IF, IS, RF, and RS, respectively. Our simulations correspond to the earliest phases of the protostars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the ideal MHD models, the properties of the protostellar cores are very similar to the spherical model. The evolution of the masses and radii of the protostellar cores are essentially identical to those in the spherical model SP. This is because the magnetic fields take almost all the angular momentum away from the gas in the central region. Contrarily, the resistive models have significantly larger angular momenta and the protostellar cores are strongly supported by rotation. To clarify the differences of the evolution between the models, we show the evolution of the radii, masses, and angular momenta of the protostellar cores in Figure 22. Here, we define the radius of the protostellar core as the radius where the infall velocity is largest (corresponding to the shock at the surface of the protostellar core) in the xy-plane, and measure the mass and angular momentum within this radius. The angular momenta in the protostellar cores in the resistive models are larger than those in the ideal MHD models by more than two orders of magnitude. The rotationally supported protostellar cores (or "circumstellar disks") are quickly built up within ∼1 year after the formation of the protostellar cores. At the end of the simulations, the radii of the disks are about 0.35 AU in RF and 0.2 AU in RS, and they are still growing as the gas with larger angular momentum accretes.

Zoom In Zoom Out Reset image size

The protostellar cores in the resistive models also look nearly spherical, but this is actually just a coincidence. They are supported by rotation in the horizontal direction, but they are vertically inflated due to the outflows. The toroidal magnetic fields are rapidly amplified in these rotating cores and the magnetic pressure gradient force drives the well-collimated outflows, or "jets" (Figures 23–25). The outflows are visible in the cross sections of the temperature as the hot towers. In the fast rotating model RF, the maximum outflow velocity reaches v_z ∼ 15 km s⁻¹, while it is ∼6 km s⁻¹ in RS. These velocities are comparable to the rotational velocities seen in the protostellar cores, and therefore far faster than the outflows driven from the first cores because of the deeper gravitational potential.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The magnetic fields in the rotating protostellar cores are quickly wound up and form the so-called magnetic tower. The tightly wound magnetic fields are susceptible to the kink instability in long wavelengths. In our resistive models, the kink instability grows rapidly and the outflows start precession (Figures 23–25). Although this instability disturbs the coherent toroidal magnetic fields, the outflow velocity is still accelerating because the bulk angular momentum in the protostellar core is increasing (Figure 22) as the matter with higher angular momentum continuously accretes from the envelope, or the remnant of the first core. Therefore, we expect that the outflow will be faster as the disk acquires larger angular momentum and the gravitational potential becomes deeper.

In all the models, radiation feedback from the protostellar core formation on the first core and outer envelope is not significant. The first cores seem to remain almost unaffected after the formation of the protostellar cores. This is different from previous RHD simulations without magnetic fields; the surrounding first cores are heated up by radiation and bipolar outflows are launched (Bate 2010, 2011; Schönke & Tscharnuter 2011). We need to study longer evolution to address this problem because we only calculate the earliest (t ≲ 1 yr) evolution of the protostellar cores; Bate (2010, 2011) and Schönke & Tscharnuter (2011) followed the evolution of the protostellar cores more than 15 years, longer than the free-fall time of the first cores.

5.3. Summary of Results

We performed 3D nested-grid RMHD simulations of formation of protostellar cores from molecular cloud cores with and without Ohmic dissipation, and revealed the earliest evolution (only one year after the formation) of the protostellar cores. These simulations are, to our knowledge, the first 3D RMHD simulations in the world following the whole evolution from molecular cloud cores to protostellar cores with realistic physics. We successfully revealed the realistic evolution in the early phase of protostellar collapse.

Previous studies (Tomida et al. 2010a, 2010b; Bate 2011) showed that the barotropic approximation fails to reproduce the realistic thermal properties, and here we also revealed that the discrepancy is more prominent in the protostellar core phase, because the barotropic approximation does not take account of the shock heating and therefore tend to underestimate the temperature, which results in the smaller radii of the first and protostellar cores and cannot reproduce the transient expansion of the protostellar cores. In order to capture the realistic thermal evolution, RMHD simulations are thus essential. Our results can be used as the initial and boundary conditions for the stellar evolution simulations (e.g., Baraffe et al. 2009; Hosokawa et al. 2011).

We found two qualitatively different phenomena between the ideal MHD models and the resistive MHD models: formation of circumstellar disks and outflows from protostellar cores. In the resistive models, the circumstellar disks are quickly built up in the vicinity of the formed protostars. Although the formed disks are still very small (R < 0.35 AU), we expect that the disks will grow continuously as the gas with larger angular momentum accretes. Machida et al. (2011) performed resistive MHD simulations with the sink cell technique and the barotropic approximation and showed that the rotationally supported disk grows to R > 200 AU while the high-density region where Ohmic dissipation works expands as the disk grows. The two-component outflows are launched from different scales as the natural by-products: the slow loosely collimated outflow driven from the first core and the fast well-collimated jet driven from the protostellar core. These are consistent with previous studies using the barotropic approximation (Machida et al. 2006a, 2007, 2008a). Although the maximum velocity of the outflow is still not very fast (v_z ∼ 15 km s⁻¹), we expect that it will get faster as the protostellar core grows. Meanwhile, we expect that the outflow from the first core is continuously driven from the outer disk or the remnant of the first core until the accretion stops (Machida & Matsumoto 2011; Machida et al. 2011). If this is the case, our simulations can naturally explain the observed high-velocity jets, especially the multi-component outflows (for example, Lee et al. 2000; Santiago-García et al. 2009). In the ideal MHD cases, on the other hand, the protostellar cores are almost the same with the spherically symmetric case and no outflow is launched as a result of efficient magnetic braking. These results do not necessarily mean that circumstellar disks are never formed in ideal MHD simulations. The circumstellar disks and outflows may be formed later when the gas with large angular momentum accretes sufficiently, because the magnetic braking is not a process that reduces the total angular momentum but transports the angular momentum from the disk to the thin envelope within the molecular cloud via Alfvén waves. To simulate such a system, long-term simulations with reasonable boundary conditions are required.

Although when and how circumstellar disks are formed is a crucial question, there has been no conclusive answer. Our results suggest that a circumstellar disk is formed in the earliest phase of protostar evolution, essentially in parallel. Jørgensen et al. (2009) suggested that Class-0 sources (young, embedded) are associated with more massive disks than Class-I sources, which support the early formation of circumstellar disks. We need more elaborate observations with higher spatial resolution to answer this question, and we expect that the Atacama Large Millimeter/submillimeter Array will provide good clues to this problem.

In our simulations, the formation of the protostellar core does not affect the first core and outer envelope. Considering the large energy released in the second collapse and subsequent accretion with the high accretion rate, we expect that the feedback like Bate (2011) and Schönke & Tscharnuter (2011) reported may happen; bipolar outflows may be launched by radiation heating from the protostellar core. There are two possible interpretations of the discrepancy: magnetic fields and duration of the simulations. Because of the efficient angular momentum transport due to magnetic fields, the first cores are not supported by rotation even in the resistive models and the structures of the accretion flows are completely different between the models with and without magnetic fields. Indeed, Bate (2011) showed that the outflows are less energetic in the models with slower rotation. Therefore, magnetic fields possibly can prevent formation of such outflows driven by radiation heating from the protostellar cores.

The other point is that we only calculate about one year after the protostellar core formation while Bate (2011) and Schönke & Tscharnuter (2011) followed the evolution longer than 15 years. Because one year is shorter than the free-fall time of the first core, the protostellar cores are still deeply obscured in the remnants of the first cores in our simulations. If this is the case, how long will it take for the feedback to become prominent? As we mentioned above, the answer is the free-fall timescale of the first core, which is about five years. Another interesting factor is the interaction between the first core and the outflow from the protostellar core. When the outflow from the protostellar core penetrates the first core, it will have some influence on the first core. Assuming the constant traveling velocity of 15 km s⁻¹, it will reach the surface of the first core in about 1.5 years. Unfortunately, it is almost hopeless to follow such a long evolution with our current code because the simulation timesteps are too small, less than one minute, and are still getting smaller. In order to perform the long-term simulations, we have to modify our code, for example, by introducing the sink particle technique (Bate et al. 1995; Krumholz et al. 2004; Federrath et al. 2010) to replace the protostellar core with a subgrid stellar evolution model (e.g., Baraffe et al. 2002; Yorke & Bodenheimer 2008; Hosokawa & Omukai 2009). To construct reasonable subgrid models, our results can be used as templates. Since the radiation feedback from the formed protostar can be highly anisotropic ("flashlight effect," Yorke & Bodenheimer 1999) due to the small-scale optically thick disk (Vaidya et al. 2009; Tanaka & Nakamoto 2011), high resolution down to the protostellar core scale is of crucial importance even when we introduce the sink particle and subgrid stellar evolution model.

All the physical processes involved in this work play crucial roles in star formation processes. Magnetic fields efficiently extract the angular momentum and drive the outflows from the first cores. Ohmic dissipation suppresses the angular momentum transport which is too efficient in the ideal MHD models and enables the formation of the circumstellar disks and the fast outflows from the protostellar cores. The EOS and radiation transfer give the realistic thermal evolution and enable us to discuss the dynamics and thermodynamics quantitatively, which cannot be achieved with the barotropic approximation. We should note, however, that there is no established test for such a complicated system, although we tested each part of our code separately with standard tests. To verify the validity of the codes, we are planning to perform comparisons with other groups.

We thank Professor Matthew R. Bate, Professor Shu-ichiro Inutsuka, Professor Kazuyuki Omukai, and Dr. Takashi Hosokawa for fruitful discussions. We also thank Professor James M. Stone for reading the manuscript carefully. We are grateful to Professor Masahiro Ikoma for providing his notes on EOS and Dr. Takahiro Miyoshi for providing information about the HLLD and HLLD− solvers. We also appreciate Richard B. Reichart checking the manuscript and improving our English. Numerical computations were partly performed on NEC SX-9 at Center for Computational Astrophysics of National Astronomical Observatory of Japan, at Japan Aerospace Exploration Agency and at Osaka University. This work is partly supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Grants-in-Aid for Scientific Research, 21244021 (Tomisaka) and 23540270 (Matsumoto). K. Tomida, Y. Hori, and S. Okuzumi are supported by the Research Fellowship from the Japan Society for the Promotion of Science (JSPS).

In some previous studies, the simple EOS of perfect gas has been often used, i.e., the adiabatic index Γ is assumed to be a constant throughout simulations, 7/5 which is valid for completely idealized diatomic molecular gas, or 5/3 which is correct only when the gas temperature is very low. However, the adiabatic index is not constant in reality. In the low temperature molecular cloud (T ≲ 100 K), molecular hydrogen behaves like monoatomic gas, Γ ∼ 5/3, because the collisional energy is insufficient to excite the rotational degrees of freedom and only translations are excited. When the gas temperature exceeds T ≳ 100 K, the rotational degrees become excited and contribute to the heat capacity. Then the adiabatic index decreases, Γ ∼ 7/5, close to that of the diatomic gas. The adiabatic index is of critical importance in the thermal evolution of the gas, and also in the stability of the gas against the gravitational instability (Stamatellos & Whitworth 2009); the stiffer gas (with larger Γ) is more stable gravitationally because it reacts more strongly against compression. Therefore, we need to consider these quantum thermodynamic behaviors in our EOS.

The second collapse is driven by the endothermic reaction of hydrogen molecule dissociation. In order to simulate the evolution in the second collapse phase, we need to take the chemical reactions into account. However, it requires quite a large computational cost to solve the network of chemical reactions in every fluid element while solving (radiation) hydrodynamics. Fortunately, because we are mainly interested in dense gas, we can assume that the timescale of chemical reactions is shorter than the dynamical timescale. Therefore, we implement the chemical reactions related to major species within the EOS on the assumption of the local thermodynamic and chemical equilibrium.

We require some assumptions to calculate the EOS for simplicity.

• 1.  
  The gas is in the local thermodynamic and chemical equilibrium (except for the ortho/para ratio of molecular hydrogen; see below).
• 2.  
  All the atoms, molecules, and ions are in the ground states.
• 3.  
  The EOS can be calculated by simple summation of each component and each degree of freedom, i.e., the interactions between components and other non-ideal effects are neglected.
• 4.  
  Only seven major species (H₂, H, H⁺, He, He⁺, He^{2 +}, and e⁻) are considered and other elements are neglected.

Based on these assumptions, we calculate the EOS using the statistical mechanics theory. Here we also assume that the gas has the solar abundance, X = 0.7 and Y = 0.28.

Partition Functions

Here, we describe partition functions of each element. The partition function of a species i can be divided into five parts; translation Z_{tr, i}, rotation Z_{rot, i}, vibration Z_{vib, i}, spin Z_{spin, i}, and electron excitation Z_{elec, i} (we include contributions from H₂ dissociation and ionization of hydrogen and helium in the electron excitation part):

$$\begin{eqnarray} &&Z_i=V\times Z_{{\rm tr},i}\times Z_{{\rm rot},i}\times Z_{{\rm vib},i} \times Z_{{\rm spin},i} \times Z_{{\rm elec},i},\quad\quad \end{eqnarray} \tag{ A1 }$$

where V is the volume. We calculate these partition functions by the standard procedure. In the following descriptions, m_i denotes the mass of i-species, h the Planck constant, and k the Boltzmann constant, respectively. The partition functions not explicitly described are unity.

Molecular hydrogen:

$$\begin{eqnarray} Z_{\rm tr,H_2}&=&\frac{(2\pi m_{\rm H_2}kT)^{3/2}}{h^3}, \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} Z_{\rm rot,H_2}&=&\left(Z_{\rm rot,H_2}^{\rm even}\right)^{{1}/{4}}\left[3 Z_{\rm rot,H_2}^{\rm odd}\exp \left(\frac{\theta _{\rm rot}}{T}\right)\right]^{{3}/{4}}, \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} Z_{\rm rot,H_2}^{\rm even}&=&\sum _{j=0,2,4,\ldots }(2j+1)\exp \left[-\frac{j(j+1)\theta _{\rm rot}}{2T}\right], \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} Z_{\rm rot,H_2}^{\rm odd}&=&\sum _{j=1,3,5,\ldots }(2j+1)\exp \left[-\frac{j(j+1)\theta _{\rm rot}}{2T}\right], \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} Z_{\rm vib,H_2}&=&\frac{1}{2\sinh {(\theta _{\rm vib}/2T)}}, \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} Z_{\rm spin,H_2}&=&\left(2\cdot \frac{1}{2}+1\right)^2=4, \end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray} Z_{\rm elec,H_2}&=&2, \end{eqnarray} \tag{ A8 }$$

where θ_rot = 170.64 K is the excitation temperature of rotation and θ_vib = 5984.48 K is that of vibration. Molecular hydrogen has two forms: orthohydrogen with aligned nuclear spins and odd rotational states, and parahydrogen with antiparallel nuclear spins and even rotational states. Here, we assumed that the ratio of orthohydrogen to parahydrogen is 3:1. This ratio has a significant impact on the dynamics of collapsing molecular cloud cores in the relatively low temperature region because thermodynamic properties related to rotation of molecular hydrogen depend on the nuclear spins and rotational states. Unfortunately, this ratio in star-forming regions is quite unclear because of observational difficulties. Although parahydrogen is more stable, some observations of interstellar dark clouds suggest that the ratio is considerably far from the equilibrium value even in the cold environment; Pagani et al. (2011) proposed that the ortho/para ratio is larger than 0.1. On the other hand, Dislaire et al. (2012) claimed that the ratio is quite small, ∼10⁻³. In this work, we calculate the EOS using the ortho/para ratio of 3:1. This assumption helps us interpret our simulation results because the adiabatic index Γ decreases monotonically (Boley et al. 2007) and also compare our results with recent simulations performed by Bate (2010, 2011; but Stamatellos & Whitworth 2009 assumed the equilibrium ratio).

Atomic hydrogen:

$$\begin{eqnarray} Z_{\rm tr,H}&=&\frac{(2\pi m_{\rm H}kT)^{3/2}}{h^3}, \end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray} Z_{\rm spin,H}&=&2\cdot \frac{1}{2}+1=2, \end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray} Z_{\rm elec,H}&=&2\exp \left(-\frac{\chi _{\rm dis}}{2kT}\right), \end{eqnarray} \tag{ A11 }$$

where χ_dis = 7.17 × 10⁻¹² erg is the dissociation energy of H₂ (Liu et al. 2009).

Ionized hydrogen:

$$\begin{eqnarray} Z_{\rm tr,H^+}&=&\frac{(2\pi m_{\rm H^+}kT)^{3/2}}{h^3}, \end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray} Z_{\rm spin,H^+}&=&2\cdot \frac{1}{2}+1=2, \end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray} Z_{\rm elec,H^+}&=&2\exp \left(-\frac{\chi _{\rm dis}+2\chi _{\rm ion}}{2kT}\right), \end{eqnarray} \tag{ A14 }$$

where χ_ion = 2.18 × 10⁻¹¹ erg is the ionization energy of atomic hydrogen.

Atomic and ionized helium:

$$\begin{eqnarray} Z_{\rm tr,He}&=&\frac{(2\pi m_{\rm He}kT)^{3/2}}{h^3}, \end{eqnarray} \tag{ A15 }$$

$$\begin{eqnarray} Z_{\rm tr,He^+}&=&\frac{(2\pi m_{\rm He^+}kT)^{3/2}}{h^3}, \end{eqnarray} \tag{ A16 }$$

$$\begin{eqnarray} Z_{\rm elec,He^+}&=&\exp \left(-\frac{\chi _{\rm He,1}}{kT}\right), \end{eqnarray} \tag{ A17 }$$

$$\begin{eqnarray} Z_{\rm tr,He^{2+}}&=&\frac{(2\pi m_{\rm He^{2+}}kT)^{3/2}}{h^3}, \end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray} Z_{\rm elec,He^{2+}}&=&\exp \left(-\frac{\chi _{\rm He,1}+\chi _{\rm He,2}}{kT}\right), \end{eqnarray} \tag{ A19 }$$

where χ_{He, 1} = 3.94 × 10⁻¹¹ erg and χ_{He, 2} = 8.72 × 10⁻¹¹ erg are the first and second ionization energies of helium.

Electron:

$$\begin{eqnarray} Z_{{\rm tr},e}&=&\frac{(2\pi m_{\rm e}kT)^{3/2}}{h^3}, \end{eqnarray} \tag{ A20 }$$

$$\begin{eqnarray} Z_{{\rm spin},e}&=&2\cdot \frac{1}{2}+1=2. \end{eqnarray} \tag{ A21 }$$

Chemical Reactions and Number Densities

The grand canonical partition function of each species is defined as

$$\begin{eqnarray} &&\Theta _i(\mu _i,V,T)=\sum _{N_i}\exp \left(\frac{N_i\mu _i}{kT}\right)\frac{Z_i^{N_i}}{N_i!}=\exp \left[Z_i\exp \left(\frac{\mu _i}{kT}\right)\right],\nonumber\\ \end{eqnarray} \tag{ A22 }$$

where μ_i is the chemical potential of i-species and N_i is the total number of i-species. The grand potential can be derived from the grand canonical partition function:

$$\begin{eqnarray} &&\hspace{-14pt}\Omega _i(\mu _i,V,T)=-kT\log \Theta _i=-kTZ_i\exp \left(\frac{\mu _i}{kT}\right). \end{eqnarray} \tag{ A23 }$$

The total grand potential can be calculated from the summation of each component:

$$\begin{eqnarray} &&\Omega =\sum _i\Omega _i. \end{eqnarray} \tag{ A24 }$$

We calculate required thermodynamic variables from these functions. First, we calculate the number density of each species based on the chemical equilibrium. The number density of i-species is derived from the partial derivative of the grand potential with respect to μ_i:

$$\begin{eqnarray} &&n_i=\frac{1}{V}\left(\frac{\partial \Omega }{\partial \mu _i}\right)_{\mu _j,V,T}=z_i\exp \left(\frac{\mu _i}{kT}\right), \end{eqnarray} \tag{ A25 }$$

where z_i = Z_i/V. This relation yields

$$\begin{eqnarray} &&\mu _i=kT\log \frac{n_i}{z_i}. \end{eqnarray} \tag{ A26 }$$

We consider only four reactions between the seven species dominant in relatively dense (but not too dense) gas for star formation problems:

$$\begin{eqnarray} {\rm H_2} & \longleftrightarrow & {\rm 2H}, \end{eqnarray} \tag{ A27 }$$

$$\begin{eqnarray} {\rm H} & \longleftrightarrow & {\rm H^+}+ {e^-}, \end{eqnarray} \tag{ A28 }$$

$$\begin{eqnarray} {\rm He} & \longleftrightarrow & {\rm He^+}+{e^-}, \end{eqnarray} \tag{ A29 }$$

$$\begin{eqnarray} {\rm He^+} & \longleftrightarrow & {\rm He^{2+}}+{e^-}. \end{eqnarray} \tag{ A30 }$$

Then the number densities can be calculated from the balance between the chemical potentials in these chemical reactions:

$$\begin{eqnarray} \mu _{\rm H_2}=2\mu _{\rm H} & \Longrightarrow & \frac{n_{\rm H}^2}{n_{\rm H_2}}=\frac{z_{\rm H}^2}{z_{\rm H_2}}\equiv K_{\rm dis}, \end{eqnarray} \tag{ A31 }$$

$$\begin{eqnarray} \mu _{\rm H}=\mu _{\rm H^+}+\mu _e & \Longrightarrow & \frac{n_{\rm H^+}n_e}{n_{\rm H}}=\frac{z_{\rm H^+}z_e}{z_{\rm H}}\equiv K_{\rm ion}, \end{eqnarray} \tag{ A32 }$$

$$\begin{eqnarray} \hspace{-12pt}\mu _{\rm He}=\mu _{\rm He^+}+\mu _e & \Longrightarrow & \frac{n_{\rm He^+}n_e}{n_{\rm He}}=\frac{z_{\rm He^+}z_e}{z_{\rm He}}\equiv K_{\rm He,1}, \end{eqnarray} \tag{ A33 }$$

$$\begin{eqnarray} \hspace{-16pt}\mu _{\rm He^+}=\mu _{\rm He^{2+}}+\mu _e & \Longrightarrow & \frac{n_{\rm He^{2+}}n_e}{n_{\rm He^+}}=\frac{z_{\rm He^{2+}}z_e}{z_{\rm He^+}}\equiv K_{\rm He,2}. \end{eqnarray} \tag{ A34 }$$

The right-hand side term of each equation, K_*, can be calculated from the partition functions. We have three additional relations; conservation of the total number density of hydrogen, conservation of the total number density of helium, and the charge neutrality:

$$\begin{eqnarray} 2n_{\rm H_2}+n_{\rm H}+n_{\rm H^+}&=&n_{\rm tot}^{\rm H}\hspace{5.0pt}\left(=\frac{\rho X}{m_{\rm H}}\right), \end{eqnarray} \tag{ A35 }$$

$$\begin{eqnarray} n_{\rm He}+n_{\rm He^+}+n_{\rm He^{2+}}&=&n_{\rm tot}^{\rm He}\hspace{5.0pt}\left(=\frac{\rho Y}{m_{\rm He}}\right), \end{eqnarray} \tag{ A36 }$$

$$\begin{eqnarray} n_{\rm H^+}+n_{\rm He^+}+2n_{\rm He^{2+}}&=&n_e. \end{eqnarray} \tag{ A37 }$$

We eliminate $n_{\rm H_2}, n_{\rm H^+}, n_{\rm He^+},$ and $n_{\rm He^{2+}}$ from Equations (A35)–(A37) using Equations (A31)–(A34):

$$\begin{eqnarray} 2\frac{n_{\rm H}^2}{K_{\rm dis}}+n_{\rm H}+\frac{n_{\rm H}}{n_e}K_{\rm ion}&=&n_{\rm tot}^{\rm H}, \end{eqnarray} \tag{ A38 }$$

$$\begin{eqnarray} n_{\rm He}\left(1+\frac{K_{\rm He,1}}{n_e}+\frac{K_{\rm He,1}K_{\rm He,2}}{n_e^2}\right)&=&n_{\rm tot}^{\rm He}, \end{eqnarray} \tag{ A39 }$$

$$\begin{eqnarray} \frac{n_{\rm H}}{n_e}K_{\rm ion}+\frac{n_{\rm He}}{n_e}K_{\rm He,1}+2\frac{n_{\rm He}}{n_e^2}K_{\rm He,1}K_{\rm He,2}&=&n_e. \end{eqnarray} \tag{ A40 }$$

By substituting n_H and n_He (we can determine the solution of (A38) uniquely since all the physical variables are positive) to (A40), we obtain one nonlinear equation of n_e:

$$\begin{eqnarray} f(n_e)&=&\frac{2n_e^2n_{\rm tot}^{\rm H}K_{\rm ion}}{\sqrt{(n_e+K_{\rm ion})^2+\frac{8}{K_{\rm dis}}n_{\rm tot}^{\rm H}n_e^2}+(n_e+K_{\rm ion})}\nonumber \\ &&+\,\frac{K_{\rm He,1}n_e+2K_{\rm He,1}K_{\rm He,2}}{n_e^2+K_{\rm He,1}n_e+K_{\rm He,1}K_{\rm He,2}}n_{\rm tot}^{\rm He}n_e^2-n_e^3=0. \nonumber\\ \end{eqnarray} \tag{ A41 }$$

Then we solve this equation numerically using the bi-section method (note that the first term of f(n_e) is already modified to avoid the round-off error). Using the obtained n_e, it is straightforward to calculate the number density of each species.

In order to derive thermodynamic variables, the temperature and density derivatives of the number densities are required. We take logarithmic differentiation of Equations (A31)–(A37), then they yield

$$\begin{eqnarray} \hspace{-4pt}&&\left(\hspace{-4.25pt} \begin{array}{@{}ccccccc@{}}1&-2&0&0&0&0&0\\ \;2n_{\rm H_2}&n_{\rm H}&n_{\rm H^+}&0&0&0&0\\ 0&-1&1&1&0&0&0\\ 0&0&n_{\rm H^+}&-n_e&0&n_{\rm He^+}&2n_{\rm He^{2+}}\\ 0&0&0&1&-1&1&0\\ 0&0&0&1&0&-1&1\\ 0&0&0&0&n_{\rm He}&n_{\rm He^+}&n_{\rm He^{2+}} \end{array} \right)\left(\begin{array}{@{}c@{}}d\ln n_{\rm H_2}\\\ms d\ln n_{\rm H}\\\ms d\ln n_{\rm H^+}\\\ms d\ln n_{\rm He}\\\ms d\ln n_{\rm He^+}\\\ms d\ln n_{\rm He^{2+}}\\\ms d\ln n_e \end{array}\right) \nonumber\\ &&\quad =\left(\begin{array}{@{}c@{}}-d\ln K_{\rm dis}\\\ms n_{\rm tot}^{\rm H}d\ln n_{\rm tot}^{\rm H}\\\ms d\ln K_{\rm ion}\\\ms 0\\\ms d\ln K_{\rm He,1}\\\ms d\ln K_{\rm He,2}\\\ms n_{\rm tot}^{\rm He}d\ln n_{\rm tot}^{\rm He} \end{array}\right). \end{eqnarray} \tag{ A42 }$$

From this matrix equation we can numerically derive the required derivatives such as (∂ln n_i/∂ln T)_ρ and (∂ln n_i/∂ln ρ)_T. For (∂ln n_i/∂ln T)_ρ the right-hand side vector becomes

$$\begin{eqnarray} &&\left(\begin{array}{@{}c@{}}\frac{d\ln z_{\rm H_2}}{d\ln T}-2\frac{d\ln z_{\rm H}}{d\ln T}\\\ms 0\\\ms \frac{d\ln z_{\rm H^+}}{d\ln T}+\frac{d\ln z_e}{d\ln T}-\frac{d\ln z_{\rm H}}{d\ln T}\\\ms 0\\\ms \frac{d\ln z_{\rm He^+}}{d\ln T}+\frac{d\ln z_e}{d\ln T}-\frac{d\ln z_{\rm He}}{d\ln T}\\\ms \frac{d\ln z_{\rm He^{2+}}}{d\ln T}+\frac{d\ln z_e}{d\ln T}-\frac{d\ln z_{\rm He^+}}{d\ln T}\\\ms 0 \end{array}\right), \end{eqnarray} \tag{ A43 }$$

and for (∂ln n_i/∂ln ρ)_T

$$\begin{eqnarray} &&\left(\begin{array}{@{}c@{}}0\\\ms n_{\rm tot}^{\rm H}\\\ms 0\\\ms 0\\\ms 0\\\ms 0\\\ms n_{\rm tot}^{\rm He} \end{array}\right). \end{eqnarray} \tag{ A44 }$$

The temperature derivatives of the partition functions can be calculated analytically.

The total number density is

$$\begin{eqnarray} &&n=\sum _i n_i, \end{eqnarray} \tag{ A45 }$$

and its derivatives are

$$\begin{eqnarray} &&\left(\frac{\partial \ln n}{\partial \ln T}\right)_\rho \equiv n_T=\sum _i \frac{n_i}{n}\left(\frac{\partial \ln n_i}{\partial \ln T}\right)_\rho, \end{eqnarray} \tag{ A46 }$$

$$\begin{eqnarray} &&\left(\frac{\partial \ln n}{\partial \ln \rho }\right)_T\equiv n_\rho =\sum _i \frac{n_i}{n}\left(\frac{\partial \ln n_i}{\partial \ln \rho }\right)_T. \end{eqnarray} \tag{ A47 }$$

Thermodynamic Variables

We use the relation valid in ideal gas:

$$\begin{eqnarray} &&P=nkT\hspace{5.0pt}\left(=\frac{\rho }{\mu m_{\rm H}}kT\right), \end{eqnarray} \tag{ A48 }$$

where μ = ρ/nm_H is the mean molecular weight. The derivatives of the pressure are

$$\begin{eqnarray} \left(\frac{\partial \ln P}{\partial \ln T}\right)_\rho \equiv P_T&=&1+n_T, \end{eqnarray} \tag{ A49 }$$

$$\begin{eqnarray} \left(\frac{\partial \ln P}{\partial \ln \rho }\right)_T\equiv P_\rho &=&n_\rho. \end{eqnarray} \tag{ A50 }$$

The specific entropy of each species is derived from the grand potential:

$$\begin{eqnarray} &&\hspace{-14pt}S_i=-\frac{1}{\rho V}\left(\frac{\partial \Omega _i}{\partial T}\right)_{V,\mu _i}=\frac{kn_i}{\rho }\left(1+\frac{d\ln z_i}{d \ln T}-\frac{\mu _i}{kT}\right),\quad \end{eqnarray} \tag{ A51 }$$

and its derivatives are

$$\begin{eqnarray} \left(\frac{\partial S_i}{\partial T}\right)_\rho &=&\frac{kn_i}{\rho T}\left[\frac{d^2\ln z_i}{d \ln T^2}+\left\lbrace 1+\left(\frac{\partial \ln n_i}{\partial \ln T}\right)_\rho \right\rbrace \frac{d\ln z_i}{d \ln T}\right]\nonumber\\ &&-\frac{\partial n_i}{\partial T}\frac{\mu _i}{\rho T}, \end{eqnarray} \tag{ A52 }$$

$$\begin{eqnarray} \left(\frac{\partial S_i}{\partial \rho }\right)_T&=&\frac{kn_i}{\rho ^2}\left[\left\lbrace \left(\frac{\partial \ln n_i}{\partial \ln \rho }\right)_T-1\right\rbrace \frac{d\ln z_i}{d\ln T}-1\right]\nonumber\\ &&+\left(\frac{n_i}{\rho }-\frac{\partial n_i}{\partial \rho }\right)\frac{\mu _i}{\rho T}. \end{eqnarray} \tag{ A53 }$$

Because the last terms in these derivatives are canceled out by taking summation of species when the chemical reactions are in equilibrium, we can omit these terms. Then the total entropy and its derivatives are defined as

$$\begin{eqnarray} S&=&\sum _i S_i, \end{eqnarray} \tag{ A54 }$$

$$\begin{eqnarray} \left(\frac{\partial S}{\partial T}\right)_\rho \equiv S_T &=&\sum _i \left(\frac{\partial S_i}{\partial T}\right)_\rho, \end{eqnarray} \tag{ A55 }$$

$$\begin{eqnarray} \left(\frac{\partial S}{\partial \rho }\right)_T\equiv S_\rho &=&\sum _i \left(\frac{\partial S_i}{\partial \rho }\right)_T. \end{eqnarray} \tag{ A56 }$$

Now we can derive thermodynamic properties we require in radiation hydrodynamic simulations.

The isothermal and adiabatic sound speeds:

$$\begin{eqnarray} c_T\equiv \sqrt{\left(\frac{\partial P}{\partial \rho }\right)_T}&=&\sqrt{\frac{P}{\rho }P_\rho }, \end{eqnarray} \tag{ A57 }$$

$$\begin{eqnarray} c_S\equiv \sqrt{\left(\frac{\partial P}{\partial \rho }\right)_S}&=&\sqrt{\frac{P}{\rho }P_\rho -\frac{P}{T}\frac{P_T S_\rho }{S_T}}. \end{eqnarray} \tag{ A58 }$$

The adiabatic index:

$$\begin{eqnarray} &&\Gamma =\left(\frac{d\ln P}{d\ln \rho }\right)_S=\frac{\rho }{P}c_S^2. \end{eqnarray} \tag{ A59 }$$

The internal energy per volume and its derivative:

$$\begin{eqnarray} e_g&=&\sum _i n_i kT\frac{d\ln z_i}{d\ln T}, \end{eqnarray} \tag{ A60 }$$

$$\begin{eqnarray} \left(\frac{\partial e_g}{\partial T}\right)_\rho &=& C_V=\rho T S_T. \end{eqnarray} \tag{ A61 }$$

Note that we do not use the relation e = C_VT which is valid only in the completely ideal cases, as Boley et al. (2007) suggested (see also Black & Bodenheimer 1975).

We tabularize these thermodynamic variables as functions of (ρ, e_g) and (ρ, T) with sufficiently high resolution (Δlog ρ = 0.05, Δlog e_g = 0.025, and Δlog T = 0.02) in ρ = 10⁻²²–10 g cm⁻³ and T = 3–10⁶ K. We use this EOS table with bi-log-linear interpolation.

Comments on EOS

For the gray radiation transfer, we use the compiled tables of the Rosseland and Planck mean opacities of Semenov et al. (2003),⁸ Ferguson et al. (2005),⁹ and the Opacity Project (OP; Seaton et al. 1994).¹⁰ For dust opacities, we adopt the composite aggregate dust model of the normal abundance from Semenov's mean opacity tables. Though Semenov's tables also contain the gas opacities, we adopt the gas opacities of Ferguson et al. (2005) and Seaton et al. (1994) because Semenov's Planck mean opacity seems to be significantly lower than other opacity tables (Ferguson et al. 2005; Seaton et al. 1994) where molecular and atomic lines dominate the opacity sources. Therefore, we connect the tables of Semenov et al. (2003) and Ferguson et al. (2005) at the temperature where all the dust components evaporate. The dust evaporation temperature (which weakly depends on the gas density, but in the typical density region, T ∼ 1400–1500 K) is given in Semenov's opacity calculation code. We use Ferguson et al. (2005) in the low temperature region (log T < 4.5) and OP in the high temperature region (4.5 < log T < 6). We tabulate these tables as functions of (ρ, T) and use them with bi-log-linear interpolation.

Unfortunately, the opacity tables do not cover the whole required region. Figure 26 shows the coverage of the opacity tables in the ρ–T plane. The typical evolution track of the central gas element and the profile in the spherically symmetric collapse are also plotted. The dust opacities of Semenov et al. (2003) cover 10⁻¹⁸ < ρ (g cm⁻³) < 10⁻⁷. OP and Ferguson et al. (2005) cover −8 < log R < 1 where R = ρ/T³₆ and T₆ = T (K)/10⁶. It is not serious that the very low density region is not covered because that region is extremely optically thin and the details of dust opacities do not matter there. We simply extrapolate the opacities by taking the nearest value at the given temperature. The high-density region is far more problematic; we do not have proper opacities for the protostellar core. However, the thermal evolution in this region is dominated by chemical reactions (dissociation and ionization) and radiation transfer is of less importance there because the gas is extremely optically thick. Therefore, we dare to extrapolate the tables in the same manner.¹¹

Zoom In Zoom Out Reset image size

We show the combined opacity tables in Figure 27 as functions of ρ and T. Note that there is quite large uncertainty in opacities, especially due to the dust models such as the structure, composition, size distribution, and so on. The thermal evolution and dynamics in star formation remain qualitatively similar even when we change the dust parameters, but the observational properties such as spectral energy distributions are directly affected by the differences between the (monochromatic) opacities.

Zoom In Zoom Out Reset image size

We calculate Ohmic resistivity considering both thermal and non-thermal particles. For non-thermal processes, we adopt the table of resistivity based on the reaction network model between dust grains and chemical species constructed by Umebayashi & Nakano (2009; see also Okuzumi 2009). The table of η_NT is given as a function of the gas density ρ, the temperature T, and the ionization rate ξ. Here, we assume the typical interstellar ionization rate due to the cosmic rays, ξ_CR ∼ 10⁻¹⁷ s⁻¹, neglecting shielding by the gas. Since the attenuation depth of the cosmic rays is about 100 g cm⁻², the ionization rate will be lower in the deep interior of the first core. In this sense, and also because we neglect ambipolar diffusion, our models correspond to the lower limit of (but still highly efficient) magnetic flux loss. Because decay of radionuclides (²⁶Al is the dominant source) considerably contributes to the ionization rate, ξ_RA ∼ 7.6 × 10⁻¹⁹ s⁻¹ (Umebayashi & Nakano 2009), the effect from neglecting the shielding is as large as about an order of magnitude at most. We should mention that there is still large uncertainty in the resistivity from the grain properties such as the structure, composition, size distribution, and so on.

In order to calculate the resistivity for our simulations, it is sufficient to estimate the contribution to the thermal ionization from the species that has low ionization energy. Potassium (K) has very low ionization energy (kT_ion ∼ 4.33 eV) and sufficiently abundant to recover good coupling between gas and magnetic fields, therefore it is the most important electron-supplying species in star formation processes. Here, we calculate the resistivity due to the thermal ionization of potassium on the assumption of thermal equilibrium using the following equation:

$$\begin{eqnarray} &&\hspace{-12pt}\eta _T =7.5\times 10^9 \exp \left(\frac{25188\,{\rm K}}{T}\right)T^{-1/4} \rho ^{1/2} \hspace{10.0pt} \, \rm {cm^2 \, s^{-1}}. \end{eqnarray} \tag{ C1 }$$

We calculate the total resistivity as follows:

$$\begin{eqnarray} &&\eta ^{-1}=\eta _T^{-1}+\eta _{NT}^{-1}, \end{eqnarray} \tag{ C2 }$$

because the resistivity is inversely proportional to the ionization degree.

We show the resistivity and magnetic Reynolds number R_m ≡ v_ffλ_J/η as a function of the gas density in Figure 28. Here, we consider the Jeans length λ_J ≡ 2πc_s(3π/32Gρ)^1/2 and the free-fall velocity v_ff ≡ (GM_J/λ_J)^1/2 (where M_J = (4π/3)λ³_Jρ) as typical length and velocity scales to estimate the magnetic Reynolds number. To draw this plot, we adopt the barotropic approximation as a typical thermal evolution to calculate the gas temperature:

$$\begin{eqnarray} &&T=\max \left[10, 10\times \left(\frac{\rho }{\rho _{{\rm crit}}}\right)^{\Gamma -1}\right] \, {\rm K}, \end{eqnarray} \tag{ C3 }$$

where ρ_crit = 2 × 10⁻¹³ g cm⁻³ is the critical density and Γ = 7/5 is the adiabatic index. The resistivity steeply decreases in ρ ≳ 10⁻⁸ g cm⁻³ because of the thermal ionization of potassium. The magnetic fields are decoupled from fluid where the magnetic Reynolds number is less than unity. Our resistivity is significantly lower than the resistivity of Kunz & Mouschovias (2009; see also Kunz & Mouschovias 2010; Dapp et al. 2012). This difference mainly comes from the ionization rates; they assume that decay of ⁴⁰K (ξ_K ∼ 2.43 × 10⁻²³ s⁻¹) is the dominant ionization source in the dense region where cosmic rays cannot penetrate and neglect the contribution from short-lived species like ²⁶Al. Our resistivity is quite similar to that used in Machida et al. (2006a).

Zoom In Zoom Out Reset image size