Does Misalignment between Magnetic Field and Angular Momentum Enhance or Suppress Circumstellar Disk Formation?

Y. Tsukamoto; S. Okuzumi; K. Iwasaki; M. N. Machida; S. Inutsuka

doi:10.3847/1538-4357/aae4dc

1. Introduction

The role of magnetic field in the gravitational collapse of a molecular cloud core has generated much discussion in the context of the formation and early evolution of circumstellar disks and of the formation of outflows and jets (e.g., Basu & Mouschovias 1994, 1995; Allen et al. 2003; Mellon & Li 2008; Inutsuka et al. 2010; Machida et al. 2011; Braiding & Wardle 2012; Dapp et al. 2012; Matsushita et al. 2017). Magnetic field transfers the angular momentum from the inner rapidly rotating region of a core to its outer slowly rotating region. The angular momentum of the inner region therefore decreases. This process is known as magnetic braking. Simulations of the gravitational collapse of typically magnetized cloud cores have shown that magnetic braking is efficient and has a strong impact on the formation and evolution of circumstellar disks (Machida et al. 2007, 2014; Price & Bate 2007; Mellon & Li 2008; Dapp & Basu 2010; Li et al. 2011; Dapp et al. 2012; Tsukamoto et al. 2015a, 2017a; Tomida et al. 2015; Lewis et al. 2015; Lewis & Bate 2017). Reviews of the protostellar collapse process and of magnetic braking can be found in Inutsuka (2012) and Tsukamoto (2016), respectively.

It has been suggested that a misalignment between the magnetic field and the angular momentum of the cloud core changes the efficiency of magnetic braking (Mouschovias 1985; Matsumoto & Tomisaka 2004; Machida et al. 2006; Hennebelle & Ciardi 2009; Joos et al. 2012; Li et al. 2013; Wurster et al. 2017; Gray et al. 2018). However, there is qualitative disagreement between previous studies. Using analytic calculations, Mouschovias (1985) suggested that magnetic braking in the case with θ = 90° (the perpendicular configuration) is much stronger than with θ = 0° (parallel configuration), where θ is the angle between the magnetic field and the angular momentum.

Matsumoto & Tomisaka (2004) conducted the first 3D simulations of the gravitational collapse of a molecular cloud core with ${\boldsymbol{B}}\nparallel {{\boldsymbol{J}}}_{\mathrm{ang}}$ (where ${{\boldsymbol{J}}}_{\mathrm{ang}}$ and ${\boldsymbol{B}}$ are the angular momentum and magnetic field of the core) and found that the angular momentum of the central region is more efficiently removed when the initial magnetic field and angular momentum of the core are in the perpendicular configuration (i.e., θ = 90°). In particular, they unequivocally showed that the perpendicular component of the angular momentum relative to the magnetic field is selectively removed from the central region during the collapse of the cloud core (Figure 12 in Matsumoto & Tomisaka 2004). Ideal MHD simulations in Wurster et al. (2017) show that the binary formation is suppressed in a core with the perpendicular configuration, although the binary is formed in the core with parallel configuration, which also implies that the magnetic braking is more efficient in the perpendicular configuration.

In contrast, Hennebelle & Ciardi (2009) and Joos et al. (2012) reported that the efficiency of magnetic braking decreases as θ increases and is minimized in the perpendicular configuration (θ = 90°). This contradicts the results of Matsumoto & Tomisaka (2004) and Li et al. (2013), who also reported that the angular momentum of the central region is greater for the perpendicular than for the parallel configuration. They suggested that the transfer of angular momentum by outflows plays a key role. All these studies adopted the ideal MHD approximation (while nonetheless imposing a small uniform resistivity to resolve the numerical difficulty in the simulation in Li et al. 2013). The physical process considered is almost the same in all these studies, but the results are contradictory nevertheless.

More recently, Masson et al. (2016) proposed, by incorporating the ohmic and ambipolar diffusions, that misalignment does not affect the magnetic braking efficiency. The degree of ionization of the cloud core is low, and the ideal MHD approximation is not valid in the high-density region ( $\rho \gtrsim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3};$ Umebayashi & Nakano 1980; Nishi et al. 1991; Nakano et al. 2002). Non-ideal effects (i.e., ohmic diffusion, Hall effect, and ambipolar diffusion) should therefore be considered. The results of Masson et al. (2016) suggest that the magnetic diffusion eliminates the effect of misalignment on the magnetic braking efficiency. Their results also imply that the difference between the magnetic braking efficiency in the parallel and perpendicular configurations observed in ideal MHD simulations originates in the high-density region $\rho \gtrsim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ because magnetic diffusion becomes dynamically important only in this region (Umebayashi & Nakano 1990; Nishi et al. 1991; Nakano et al. 2002; Tomida et al. 2015; Tsukamoto et al. 2015b). We should note, however, that Masson et al. (2016) investigated only the core with θ < 40°. It is therefore unclear whether their claim also holds for θ > 40°. Wurster et al. (2017) reported that the separation of the binary formed in the core with the perpendicular configuration is less than with the parallel configuration. This implies that the magnetic braking is more efficient with the perpendicular configuration even with the non-ideal effects.

In summary, the earlier studies have reported apparently contradicting results, leaving the question open as to whether magnetic braking is enhanced by, suppressed by, or independent of misalignment.

Resolving this discrepancy is of interest not only for theorists but also for observational astronomers. Several observations of dust polarized emission have identified the direction of the magnetic field on various scales in molecular clouds (Girart et al. 2006; Hull et al. 2013, 2014, 2017a, 2017b; Davidson et al. 2014; Segura-Cox et al. 2015; Ward-Thompson et al. 2017; Alves et al. 2018; Cox et al. 2018; Galametz et al. 2018; Kwon et al. 2018a; Kwon et al. 2018b; Maury et al. 2018; Sadavoy et al. 2018; Soam et al. 2018).

Furthermore, recent high-resolution observations of Class 0/I YSOs have reported the existence of Keplerian disks with sizes of several tens to hundreds of au (Takakuwa et al. 2012; Tobin et al. 2012; Murillo et al. 2013; Chou et al. 2014; Ohashi et al. 2014; Aso et al. 2015, 2017; Yen et al. 2015, 2017). This information allows a discussion of the relative orientations of the magnetic field and the circumstellar disks and of the correlation between the disk size and magnetic field direction. Understanding the nature of magnetic braking in a misaligned system is necessary to interpret such observational results appropriately and to discuss the importance of a magnetic field on the evolution of angular momentum.

This paper aims to resolve the apparent discrepancies between the earlier theoretical studies. In particular, we focus on the dependence of magnetic braking on the gravitational stability of the initial cloud core, i.e., the degree of departure from gravitational equilibrium and the mass accretion rate onto the disk. Following convention, the gravitational stability of the cloud core is parameterized by $\alpha \equiv {E}_{\mathrm{therm}}/{E}_{\mathrm{grav}}$ , where E_therm and E_grav are the thermal and gravitational energies of the initial core, respectively. Machida et al. (2016) demonstrated that the magnetic braking efficiency increases as α increases. However, Machida et al. (2016) only investigated the core with the parallel configuration. On the other hand, previous studies concerning misalignment did not carefully consider the dependence of the result on α. We also conducted simulations with ohmic and ambipolar diffusions, to compare with the results by Masson et al. (2016).

The present paper is organized as follows. We first describe the numerical methods and initial conditions used in this study in Section 2. The simulation results are then shown in Section 3. We summarize our numerical results and discuss their implications in Section 4. Finally, we present our conclusions on the impact of misalignment in Section 5. In the Appendix, we review the analytic discussion of the magnetic braking timescale. The derived timescales are used for discussions.

2. Numerical Method and Initial Conditions

2.1. Numerical Method

The simulations solve non-ideal radiation magnetohydrodynamics equations with self-gravity:

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

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

$\begin{eqnarray}&&\displaystyle \frac{D}{{Dt}}\left(\displaystyle \frac{{E}_{{\rm{r}}}}{\rho }\right)=-\displaystyle \frac{{\rm{\nabla }}\cdot {{\boldsymbol{F}}}_{{\rm{r}}}}{\rho }-\displaystyle \frac{{\rm{\nabla }}{\boldsymbol{v}}:{{\mathbb{P}}}_{{\rm{r}}}}{\rho }+{\kappa }_{{\rm{P}}}c({a}_{{\rm{r}}}{T}_{{\rm{g}}}^{4}-{E}_{{\rm{r}}}),\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{D}{{Dt}}\left(\displaystyle \frac{e}{\rho }\right) & = & -\displaystyle \frac{1}{\rho }{\rm{\nabla }}\cdot \left\{\left(P,+,\displaystyle \frac{1}{2},| ,{\boldsymbol{B}},{| }^{2}\right){\boldsymbol{v}}-{\boldsymbol{B}}\left({\boldsymbol{B}},\cdot ,{\boldsymbol{v}}\right)\right\}\\ & & -\,{\kappa }_{{\rm{P}}}c({a}_{{\rm{r}}}{T}_{{\rm{g}}}^{4}-{E}_{{\rm{r}}})-{\boldsymbol{v}}\cdot {\rm{\nabla }}{\rm{\Phi }}\\ & & -\,\displaystyle \frac{1}{\rho }{\rm{\nabla }}\cdot [\{({\eta }_{{\rm{O}}}({\rm{\nabla }}\times {\boldsymbol{B}})\\ & & -\,{\eta }_{{\rm{A}}}(({\rm{\nabla }}\times {\boldsymbol{B}})\times \hat{{\boldsymbol{B}}})\times \hat{{\boldsymbol{B}}}\}\times {\boldsymbol{B}}],\end{array}\end{eqnarray} \tag{ 4 }$

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

where ρ is the gas density, P is the gas pressure, ${\boldsymbol{B}}$ is the magnetic field and $\hat{{\boldsymbol{B}}}\equiv {\boldsymbol{B}}/| {\boldsymbol{B}}|$ , ${\eta }_{{\rm{O}}}$ and ${\eta }_{{\rm{A}}}$ are, respectively, the resistivity for ohmic and ambipolar diffusion, E_r is the radiation energy, ${{\boldsymbol{F}}}_{{\rm{r}}}$ is the radiation flux, ${{\mathbb{P}}}_{{\rm{r}}}$ is the radiation pressure, T_g is the gas temperature, ${\kappa }_{{\rm{P}}}$ is the Plank mean opacity, $e=\rho u+\tfrac{1}{2}(\rho {{\boldsymbol{v}}}^{2}+{{\boldsymbol{B}}}^{2})$ is the total energy (where u is specific internal energy), Φ is the gravitational potential, and a_r and G are the radiation and gravitational constants, respectively. The Hall effect is not considered in this study but is discussed in Krasnopolsky et al. (2011), Li et al. (2011), Tsukamoto et al. (2015a), Wurster et al. (2016), Tsukamoto et al. (2017a), and Wurster et al. (2018a).

To close the equations for radiation transfer, we employed the flux-limited diffusion approximation (Levermore & Pomraning 1981),

$\begin{eqnarray*}\begin{array}{rcl}{{\boldsymbol{F}}}_{{\rm{r}}} & = & \displaystyle \frac{c\lambda }{{\kappa }_{{\rm{R}}}\rho }{\rm{\nabla }}{E}_{{\rm{r}}},\,\lambda (R)=\displaystyle \frac{2+R}{6+2R+{R}^{2}},\\ R & = & \displaystyle \frac{| {\rm{\nabla }}{E}_{{\rm{r}}}| }{{\kappa }_{{\rm{R}}}\rho {E}_{{\rm{r}}}},\,{{\mathbb{P}}}_{{\rm{r}}}={\mathbb{D}}{E}_{{\rm{r}}},\\ {\mathbb{D}} & = & \displaystyle \frac{1-\chi }{2}{\mathbb{I}}+\displaystyle \frac{3\chi -1}{2}{\boldsymbol{n}}\otimes {\boldsymbol{n}},\,\chi =\lambda +{\lambda }^{2}{R}^{2},\,\\ {\boldsymbol{n}} & = & \displaystyle \frac{{\rm{\nabla }}{E}_{{\rm{r}}}}{| {\rm{\nabla }}{E}_{{\rm{r}}}| },\end{array}\end{eqnarray*}$

where κ_R is the Rosseland mean opacity.

We use the smoothed particle hydrodynamics (SPH) method (Gingold & Monaghan 1977; Lucy 1977; Monaghan & Lattanzio 1985); the numerical code was developed by the present authors and employed previously (e.g., Tsukamoto & Machida 2011, 2013; Tsukamoto et al. 2013b, 2015a; Yoneda et al. 2016). The ideal MHD part is solved with the Godunov smoothed particle magnetohydrodynamics (GSPMHD) method (Iwasaki & Inutsuka 2011). The divergence-free condition is maintained with the hyperbolic divergence cleaning method for GSPMHD (Iwasaki & Inutsuka 2013). Radiative transfer is solved implicitly using the method of Whitehouse & Bate (2004) and Whitehouse et al. (2005). We treat ohmic and ambipolar diffusions with the methods of Tsukamoto et al. (2013a) and Wurster et al. (2014), respectively. Both diffusion processes were accelerated by super-time-stepping (Alexiades et al. 1996). To calculate the self-gravity, we adopted the Barnes–Hut tree algorithm with an opening angle of θ_gravity = 0.5 (Barnes & Hut 1986). For the dust and gas opacities, the data tables from Semenov et al. (2003) and from Ferguson et al. (2005), respectively, were adopted. We adopted the tabulated equation of state (EOS) used in Tomida et al. (2013), in which the internal degrees of freedom and chemical reactions of seven species ${{\rm{H}}}_{2},\,{\rm{H}},\,{{\rm{H}}}^{+},\,\mathrm{He},\,{\mathrm{He}}^{+},{\mathrm{He}}^{++},{e}^{-}$ are considered.

We used tabulated resistivity calculated by the methods described in Nakano et al. (2002) and Okuzumi (2009), in which we considered ion species of ${{\rm{H}}}_{3}^{+},\,{\mathrm{HCO}}^{+},\,{\mathrm{Mg}}^{+},{\mathrm{He}}^{+},{{\rm{C}}}^{+},{{\rm{H}}}^{+},{e}^{-}$ . The reaction rate was taken from the UMIST2012 database (McElroy et al. 2013). We took into account the cosmic-ray and thermal ionization, gas-phase and dust-surface recombination, and ion-neutral reaction. The dust size and density are a = 3.5 × 10⁻² μm and ${\rho }_{{\rm{d}}}=2\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , respectively. We assumed a dust-to-gas mass ratio of 1.7% to mimic the gas of the minimum-mass solar nebula model (Hayashi 1981). Our calculations considered nonthermal ionization by the cosmic-ray and thermal ionization of potassium. The cosmic-ray ionization rate was fixed to be ${\xi }_{\mathrm{CR}}={10}^{-17}\,{{\rm{s}}}^{-1}$ .

2.2. Initial and Boundary Conditions

We modeled the initial cloud cores as isothermal uniform gas spheres. The mass and temperature of the initial core are 1 M_⊙ and 10 K, respectively. The rotation energy normalized by the gravitational energy of the initial core is E_rot/E_grav = 0.03. The initial angular momentum vector is parallel to the z-axis. The initial mass-to-flux ratio relative to the critical value is $\mu ={(M/{\rm{\Phi }})/(M/{\rm{\Phi }})}_{\mathrm{crit}}=4$ , where ${\rm{\Phi }}=\pi {R}_{\mathrm{core}}^{2}{B}_{0}$ (R_core is the core radius and B₀ is the initial magnetic field strength) and (M/Φ)_crit = (0.53/3π)(5/G)^1/2 (Mouschovias & Spitzer 1976). The initial magnetic field is uniform and tilted relative to the z-axis and is given by

$\begin{eqnarray}&&{\boldsymbol{B}}=({B}_{x},{B}_{y},{B}_{z})={B}_{0}(-\sin \theta ,0,\cos \theta ).\end{eqnarray} \tag{ 6 }$

In the following simulations, we scanned over the parameter space comprising two parameters: the ratio α of the thermal energy to the gravitational energy ( $\alpha \equiv {E}_{\mathrm{therm}}/{E}_{\mathrm{grav}}$ ), and the relative angle θ between the magnetic field and angular momentum of the cloud core. In practice, different values of α are realized by changing the core radius, while the core mass is fixed to $1\,{M}_{\odot }$ . Thus, the rotation energy is also different between cores with different values of α. Table 1 lists the model names, α, the initial radius, the initial density, the initial angular velocity, the initial magnetic field strength, and θ for the cloud core. For comparison, we also simulated the cores with very weak magnetic field where the mass-to-flux ratio was μ = 10⁴ for α = 0.2, 0.4, and 0.6 and θ = 0°, which essentially correspond to the hydrodynamics simulations. Thus, we refer to these simulations as hydrodynamics simulations (results are shown in Figure 11). The initial cores were modeled with $3\times {10}^{6}$ SPH particles. We conducted the simulations until the epoch just after protostar formation (when the central density ρ_c becomes $\sim {10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ).

Table 1. Model Names and Simulation Parameters: α = E_therm/E_grav, the Initial Radius R₀, the Initial Density ρ₀, the Initial Angular Velocity Ω₀, the Initial Magnetic Field Strength ${{\boldsymbol{B}}}_{0}$ , the Angle between the Initial Magnetic Field and the Initial Angular Momentum Vector of Cloud Cores θ

Model Name	α	R₀ (au)	Density ${\rho }_{0}\,({\rm{g}}\,{\mathrm{cm}}^{-3})$	${{\rm{\Omega }}}_{0}\,({{\rm{s}}}^{-1})$	${{\boldsymbol{B}}}_{0}\,({\rm{G}})$	θ (deg)	Magnetic Diffusion
I1	0.2	2.0 × 10³	1.7 × 10⁻¹⁷	6.5 × 10⁻¹³	3.5 × 10⁻⁴	0°	No
I2	0.2	2.0 × 10³	1.7 × 10⁻¹⁷	6.5 × 10⁻¹³	3.5 × 10⁻⁴	45°	No
I3	0.2	2.0 × 10³	1.7 × 10⁻¹⁷	6.5 × 10⁻¹³	3.5 × 10⁻⁴	90°	No
I4	0.4	4.1 × 10³	2.0 × 10⁻¹⁸	2.3 × 10⁻¹³	8.6 × 10⁻⁵	0°	No
I5	0.4	4.1 × 10³	2.0 × 10⁻¹⁸	2.3 × 10⁻¹³	8.6 × 10⁻⁵	45°	No
I6	0.4	4.1 × 10³	2.0 × 10⁻¹⁸	2.3 × 10⁻¹³	8.6 × 10⁻⁵	90°	No
I7	0.6	6.1 × 10³	6.1 × 10⁻¹⁹	1.2 × 10⁻¹³	3.8 × 10⁻⁵	0°	No
I8	0.6	6.1 × 10³	6.1 × 10⁻¹⁹	1.2 × 10⁻¹³	3.8 × 10⁻⁵	45°	No
I9	0.6	6.1 × 10³	6.1 × 10⁻¹⁹	1.2 × 10⁻¹³	3.8 × 10⁻⁵	90°	No

NI1	0.2	2.0 × 10³	1.7 × 10⁻¹⁷	6.5 × 10⁻¹³	3.5 × 10⁻⁴	0°	Yes
NI2	0.2	2.0 × 10³	1.7 × 10⁻¹⁷	6.5 × 10⁻¹³	3.5 × 10⁻⁴	45°	Yes
NI3	0.2	2.0 × 10³	1.7 × 10⁻¹⁷	6.5 × 10⁻¹³	3.5 × 10⁻⁴	90°	Yes
NI4	0.4	4.1 × 10³	2.0 × 10⁻¹⁸	2.3 × 10⁻¹³	8.6 × 10⁻⁵	0°	Yes
NI5	0.4	4.1 × 10³	2.0 × 10⁻¹⁸	2.3 × 10⁻¹³	8.6 × 10⁻⁵	45°	Yes
NI6	0.4	4.1 × 10³	2.0 × 10⁻¹⁸	2.3 × 10⁻¹³	8.6 × 10⁻⁵	90°	Yes
NI7	0.6	6.1 × 10³	6.1 × 10⁻¹⁹	1.2 × 10⁻¹³	3.8 × 10⁻⁵	0°	Yes
NI8	0.6	6.1 × 10³	6.1 × 10⁻¹⁹	1.2 × 10⁻¹³	3.8 × 10⁻⁵	45°	Yes
NI9	0.6	6.1 × 10³	6.1 × 10⁻¹⁹	1.2 × 10⁻¹³	3.8 × 10⁻⁵	90°	Yes

Note. The last column indicates whether the ohmic and ambipolar diffusions are included ("yes") or not ("no").

Download table as: ASCII Typeset image

A boundary condition was imposed at R_out = 0.95R_core, and the particles with r > R_out rotated with the initial angular velocity. In addition, another boundary condition for radiative transfer was introduced by fixing the gas and radiation temperatures to 10 K in $\rho \lt 4.0\times {10}^{-17}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ .

3. Results

3.1. Density and Angular Momentum Structures

3.1.1. Ideal MHD Simulations

Figure 1 shows the density cross section in the y–z plane, obtained from the ideal MHD simulations. In the case of a nonzero θ, the system loses its axial symmetry about the z-axis, and the structures on the x–z and y–z planes, for example, are different. In particular, the outflow is not found to lie within the x–z plane in our results for $\theta \ne 0^\circ$ . By examining the cross section around the z-axis with 0° < ϕ < 180°, where ϕ is the azimuthal angle of the intersection of the cross section on the x–y plane and ϕ = 0° corresponds to the cross section of the x–z plane, we find that the outflow sits on or close to the y–z plane (or a plane with ϕ = 90°) with our initial conditions. Thus, we chose the y–z plane to investigate the density and velocity structures.

**Figure 1.** Density cross sections in the y–z plane for the central 600 au square region at the end epoch (at the epoch when the central density has reached ${\rho }_{{\rm{c}}}\sim {10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) in the simulations of the ideal MHD models (models I1–I9). The black and red lines represent the contours of the density and of ${v}_{{\rm{r}}}\equiv ({{yv}}_{y}+{{zv}}_{z})/\sqrt{{y}^{2}+{z}^{2}}=0$ , respectively. Black contour levels are $\rho ={10}^{-17},{10}^{-16},\,{10}^{-15},\,{10}^{-14},{10}^{-13},$ and ${10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ . The region enclosed by the red contours is the outflow. The red arrows represent the velocity field and the white arrows the direction of the magnetic field.
Download figure:
Standard image High-resolution image

Figure 1 clearly shows two notable structures: the pseudo-disk and the outflow. The pseudo-disk has the morphology of a flattened disk. The inclination of its normal vector relative to the z-axis corresponds roughly to the value of θ in the simulations. At the beginning of the isothermal collapse phase, the collapse is spherically symmetric and the central magnetic field is amplified according to B ∝ ρ^2/3. The spherical collapse continues for as long as the Lorentz force is negligible. Once the Lorentz force has become comparable to the gravity force, it begins to deflect the gas motion in the direction parallel to the magnetic field, and the spherical symmetry of the isothermal collapse breaks down. As a result, the fluid elements at a certain direction selectively accrete onto the midplane, and the flattened disk-like structure, i.e., pseudo-disk, is formed. Thus, the existence of a pseudo-disk or flattened envelope indicates that nonspherical collapse has occurred in the isothermal collapse phase (see Tsukamoto et al. 2015b, for detailed evolution of the magnetic field strength). As shown below, the nonspherical collapse affects the evolution of the angular momentum of the central region and the early evolution of the disk.

Figure 1 also shows contours of ${v}_{{\rm{r}}}\equiv ({{yv}}_{y}+{{zv}}_{z})/\sqrt{{y}^{2}+{z}^{2}}=0$ (in red). The enclosed regions are outflows, which are observed to be formed in all the simulations with θ = 0° and 45°. The outflow size decreases as α decreases because the protostar quickly forms in the simulations with a small α and there is not enough time for the outflow to grow. The outflow is generated by the rotation of the central region and has the angular momentum extracted from the central region (Tomisaka 2002; Machida et al. 2008). Figure 2 shows the spatial map of the angular momentum per unit volume ${J}_{{yz}}={({J}_{y}^{2}+{J}_{z}^{2})}^{1/2}$ . J_yz does not include the contribution of the rotation in the y–z plane. The figure shows that the outflow has greater angular momentum than the infalling region surrounding it (as is particularly clear in panels (c) and (f)). Hence, this demonstrates that the angular momentum has been extracted from the central region. The outflow size decreases as α decreases, whereas the angular momentum per unit volume increases (note that the color scale increases as α decreases). Interestingly, the outflow is found to form in model I9 (panel (i) of Figure 1 and 2). This is consistent with the spiral-flow-type outflow reported by Matsumoto et al. (2017).

**Figure 2.** Cross sections of the angular momenta per unit volume ${J}_{{yz}}={({J}_{y}^{2}+{J}_{z}^{2})}^{1/2}$ in the y–z plane for the central 600 au square region at the end epoch (at the epoch when the central density has reached ${\rho }_{{\rm{c}}}\sim {10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) in the simulations of the ideal MHD models (models I1–I9). The black and red lines represent contours of J_yz of ${v}_{{\rm{r}}}\equiv ({{yv}}_{y}+{{zv}}_{z})/\sqrt{{y}^{2}+{z}^{2}}=0$ , respectively. The region enclosed by the red contours is outflow. The black contour levels are ${J}_{{yz}}={10}^{4},{10}^{5}$ , and ${10}^{6}\,{\rm{g}}\,{\mathrm{cm}}^{-1}\,{{\rm{s}}}^{-1}$ . The red arrows show the velocity field.
Download figure:
Standard image High-resolution image

**Figure 2.** Cross sections of the angular momenta per unit volume ${J}_{{yz}}={({J}_{y}^{2}+{J}_{z}^{2})}^{1/2}$ in the y–z plane for the central 600 au square region at the end epoch (at the epoch when the central density has reached ${\rho }_{{\rm{c}}}\sim {10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) in the simulations of the ideal MHD models (models I1–I9). The black and red lines represent contours of J_yz of ${v}_{{\rm{r}}}\equiv ({{yv}}_{y}+{{zv}}_{z})/\sqrt{{y}^{2}+{z}^{2}}=0$ , respectively. The region enclosed by the red contours is outflow. The black contour levels are ${J}_{{yz}}={10}^{4},{10}^{5}$ , and ${10}^{6}\,{\rm{g}}\,{\mathrm{cm}}^{-1}\,{{\rm{s}}}^{-1}$ . The red arrows show the velocity field.
Download figure:
Standard image High-resolution image

In the cores with θ = 45° (panels (d)–(f)), the central structure is warped and its normal is not parallel to the z-axis (the direction of the initial angular momentum). This indicates that magnetic braking changes the direction of the angular momentum. Matsumoto & Tomisaka (2004) reported that the magnetic field changes the direction of the angular momentum during gravitational collapse. The formation of warped structure is consistent with their result.

3.1.2. Non-ideal MHD Simulations

Figures 3 and 4 show the density and angular momentum cross section in the y–z plane, respectively, as obtained with the non-ideal MHD simulations. The pseudo-disks are also formed in the non-ideal MHD simulations, indicating that the radial magnetic tension causes a deviation of the gas motion from the spherical symmetric collapse also in the non-ideal MHD simulations. This is well explained by the fact that magnetic diffusions do not play a role at low densities $\rho \lesssim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ (Umebayashi & Nakano 1990; Nakano et al. 2002; Tomida et al. 2015; Tsukamoto et al. 2015b).

**Figure 3.** Density cross sections in the y–z plane for the central 1000 au square region at the end epoch (at the epoch when the central density has reached ${\rho }_{{\rm{c}}}\sim {10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) in the simulations of the non-ideal MHD models (models NI1–NI9). The black and red lines represent contours of the density and of ${v}_{{\rm{r}}}\equiv ({{yv}}_{y}+{{zv}}_{z})/\sqrt{{y}^{2}+{z}^{2}}=0$ , respectively. The black contour levels are $\rho ={10}^{-17},{10}^{-16},\,{10}^{-15},\,{10}^{-14},{10}^{-13},$ and ${10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ . The region enclosed by the red contours is outflow. The red arrows represent the velocity field and the white arrows the direction of the magnetic field.
Download figure:
Standard image High-resolution image

**Figure 4.** Cross section of the angular momenta per unit volume ${J}_{{yz}}={({J}_{y}^{2}+{J}_{z}^{2})}^{1/2}$ in the y–z plane for the central 1000 au square region at the end epoch (at the epoch when the central density has reached ${\rho }_{{\rm{c}}}\sim {10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) of the non-ideal MHD simulations (models NI1–NI9). The black and red lines represent contours of the density and ${v}_{{\rm{r}}}\equiv ({{yv}}_{y}+{{zv}}_{z})/\sqrt{{y}^{2}+{z}^{2}}=0$ , respectively. The region enclosed by the red contours is the outflow. The black contour levels are ${J}_{{yz}}={10}^{4},{10}^{5}$ , and ${10}^{6}\,{\rm{g}}\,{\mathrm{cm}}^{-1}\,{{\rm{s}}}^{-1}$ . The red arrows represent the velocity field.
Download figure:
Standard image High-resolution image

An outflow appears in the simulations with θ = 0° and 45°. The inclusion of the magnetic diffusions in the simulation results in the weaker outflow activity. The outflow velocity and angular momentum per unit volume (Figure 4) are smaller than those in the ideal MHD simulations. Their size is, however, larger than in the ideal MHD simulations. This reflects the fact that, in the non-ideal MHD simulations, the time span from the first-core formation up to the second collapse is longer than in the ideal MHD simulations, owing to the rotation support. Thus, the system has a longer amount of time to develop a larger outflow (see Figure 7).

Figures 2 and 4 indicate that the morphology of the angular momentum map is similar for the ideal and non-ideal MHD simulations, although strengths differ. On the other hand, the outflow morphology is different. For example, panels (d) and (e) of Figures 2 and Figure 4 show that the gas in the central region of $| y| ,| z| \lt 100\,\mathrm{au}$ in non-ideal MHD simulations is not outflowing, whereas that of the ideal MHD simulation is outflowing. Panels (i) of Figures 2 and 4 show that a spiral-flow-type outflow does not appear in the non-ideal MHD simulations. These results show that the outflow activity depends on the strength of the rotation in the envelope and on whether or not magnetic diffusion is considered (see also Tsukamoto et al. 2015b; Masson et al. 2016, for the relation between the outflow activity and magnetic diffusions).

3.2. Evolution of the Angular Momentum of the Central Region

3.2.1. Ideal MHD Simulations

Figure 5 shows the evolution of the mean angular momentum of the central region (the region with $\rho \gt {10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) of ideal MHD simulations. The horizontal axis represents the central density ρ_c. Because ρ_c is an increasing function of time (see Figure 7), Figure 5 shows the time evolution of the mean angular momentum around the center.

In the models with α = 0.6 (panel (c)), which is a similar value to that adopted in Matsumoto & Tomisaka (2004), the angular momentum at the end of the simulation (at the epoch when the central density has reached ${\rho }_{{\rm{c}}}\sim {10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) is a decreasing function of the relative angle θ. This suggests that magnetic braking is stronger in a cloud core with θ = 90° than θ = 0°. This result is consistent with Matsumoto & Tomisaka (2004) and with the analytic argument of Mouschovias (1985).

As α decreases, the difference in the angular momenta becomes smaller, and in the models with α = 0.2, the central angular momentum of the core with θ = 90° is slightly greater than θ = 0° (right panel). This is consistent with the results of Hennebelle & Ciardi (2009) and Joos et al. (2012), in which the initial cloud core had α = 0.25. We note that, as in Joos et al. (2012), the difference of the central angular momentum is small at the epoch of the protostar formation. This suggests that the evolution of the central angular momentum in the misaligned cloud core depends on α when the ideal MHD approximation is adopted.

The figure shows that the mean angular momentum increases in the density range of ${10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}\lesssim {\rho }_{{\rm{c}}}\lesssim {10}^{-10}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ . This increase is caused by the mass and angular momentum accretion from the envelope. Then, it begins to decrease in the range of ${10}^{-10}\,{\rm{g}}\,{\mathrm{cm}}^{-3}\lesssim {\rho }_{{\rm{c}}}\lesssim {10}^{-8}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ . The decrease of the central angular momentum is caused by the removal of angular momentum from the central region as a result of magnetic braking. This is significant in the simulations with a large α because the mass and angular momentum accretion rates decrease as α increases and the removal of angular momentum by the magnetic field overtakes the angular momentum supply in the simulations with a greater α value.

By comparing the simulations with the same θ values (i.e., lines of the same colors in the panels of Figure 5), we observe that the central angular momentum is a decreasing function of α. This result suggests that the circumstellar disk forms more easily in the core when α is small, consistent with previous studies (Machida et al. 2014, 2016).

3.2.2. Non-ideal MHD Simulations

Figure 6 shows the evolution of the mean angular momentum of the central region with $\rho \gt {10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ in the non-ideal MHD simulations. In all the models, the central angular momentum decreases as θ increases. In particular, comparing the results between the cases of θ = 0° and θ = 90°, we find that the central angular momentum with θ = 0° is greater than with θ = 90°, independently of α.

Because the magnetic diffusions are effective in the relatively high density of $\rho \gtrsim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , the suppression of magnetic braking in the high-density region causes the difference between the ideal and non-ideal MHD simulations. The decrease of the angular momentum for ${\rho }_{{\rm{c}}}\gt {10}^{-10}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ as observed in the ideal MHD simulations (Figure 5) does not appear once magnetic diffusions have been considered. This result indicates that the removal of angular momentum from the central region is not significant and that supply of angular momentum by mass accretion dominates its removal.

Figure 6 also shows that the differences in the central angular momentum between the models with θ = 0° and θ = 45° are small. For all values of α, the central angular momenta differ between the models with θ = 0° and θ = 45° by factor of two or less. This is also true in the ideal MHD simulations (see Figure 5). Thus, we conclude that a misalignment with a small θ (θ < 45°) does not significantly affect the evolution of the central angular momentum regardless of whether non-ideal effects are included or not. The small difference in angular momentum with θ < 45° is consistent with a previous study (Masson et al. 2016), where the authors investigated the evolution of the angular momentum in the core with θ < 40° and reported that the misalignment does not introduce any significant difference in the angular momentum in the disk if magnetic diffusions are taken into account.

In the non-ideal MHD simulations, the central angular momentum is also a decreasing function of α for a given θ, consistent with the results of ideal MHD simulations. This suggests that the difference in α mainly affects magnetic braking in the isothermal collapse phase ( ${\rho }_{{\rm{c}}}\lt {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ), in which the magnetic diffusion does not play a role.

3.2.3. Evolution of the Central Density

In Figures 5 and 6, the central density was chosen as a variable describing temporal evolution. A similar parameterization has been commonly used (e.g., Matsumoto & Tomisaka 2004; Bate et al. 2014; Tomida et al. 2015). This allows us to compare simulations at a similar evolutionary stage, for example, at the epochs of the first- and ( ${\rho }_{{\rm{c}}}\sim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) and second-core (or protostar) formation ( ${\rho }_{{\rm{c}}}\sim {10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ). However, the elapsed time of each evolutionary epoch differs between the simulations. We here investigate the difference in the temporal evolution of the central density.

Figure 7 shows the temporal evolution of the central density. The time at the second collapse is greater than the freefall time by a factor of between 1.1 and 1.4. where the freefall time is t_ff = 1.6 × 10⁴ yr for α = 0.2, t_ff = 4.6 × 10⁴ yr for α = 0.4, and t_ff = 8.5 × 10⁴ yr for α = 0.6. The figure shows that the epochs of the first- ( ${\rho }_{{\rm{c}}}\sim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) and second-core formations ( ${\rho }_{{\rm{c}}}\sim {10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) depend on θ, even for the same α. In the ideal MHD simulations, the models with a greater θ require a longer time for the first- and second-core formations. For example, in the ideal MHD simulations with α = 0.6, the second-core formation epoch differs by ∼5 × 10³ yr between the simulations with θ = 0° and θ = 90°. This introduces complexity when we compare the simulations that have different θ values even when their initial cores had the same α and β_rot.

Figure 7 also shows the evolution of the central density in the non-ideal MHD simulations. For all the values of α, the epoch of the first-core formation ( ${\rho }_{{\rm{c}}}\sim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) is delayed as θ increases, consistent with the results of the ideal MHD simulations. Note that the temporal evolution of ρ_c is almost identical in the isothermal collapse phase ${\rho }_{{\rm{c}}}\lesssim {10}^{-14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ between the ideal and non-ideal MHD simulations.

The epoch of the second-core formation in the non-ideal MHD simulations increases as θ increases, for α = 0.4 and 0.6, although the difference between the result with θ = 0° and θ = 90° becomes smaller than that of the ideal MHD simulations. In the simulation with α = 0.2, it is longer for θ = 0° than for θ = 90°. The reason for the small difference of the second-collapse epoch when α = 0.4 and 0.6, and for why the second collapse happens earlier in the simulation with θ = 90° when α = 0.2, is the rotational support. Figure 6 shows that the central region for the models with θ = 0° and 45° has a much greater angular momentum than with θ = 90°. The centrifugal radius ${r}_{\mathrm{cent}}$ is estimated to be ${r}_{\mathrm{cent}}\sim {j}^{2}/{({GM})\sim 4.5(j/(3\times {10}^{19}{\mathrm{cm}}^{2}{{\rm{s}}}^{-1}))}^{2}(M/(0.1\,{M}_{\odot }))$ au. The centrifugal force with a specific angular momentum of the order of $j\sim {10}^{19}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ provides the rotational support against self-gravity at a radius of several au, which corresponds to the radius of the first core. Figures 6 and 7 show that the duration of the first-core epoch ( ${10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}\lt {\rho }_{{\rm{c}}}\,\lt {10}^{-8}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) is considerably longer in the models with $j\gtrsim {10}^{19}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ than those with $j\lesssim {10}^{19}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ . The difference in the second-collapse epoch arises in consequence.

It is worth noting that the order reversal of the second-collapse epoch observed in our simulations has a different physical origin from that observed in Vaytet et al. (2018), in which the second collapse in the ideal MHD simulation happened later than in the non-ideal MHD simulation. As those authors discussed, the order reversal they observed was caused by interchange instability.

3.2.4. Summary of Section 3.2

This subsection demonstrated that the evolution in angular momentum of the central region depends on the α value of the initial cloud cores and on whether magnetic diffusions are taken into account. In particular, it was seen that even the qualitative dependence of the central angular momentum on θ varies, depending on α in ideal MHD simulations. For a small value of α such as α = 0.2, the central angular momentum of the simulation with θ = 0 is smaller than that with θ = 90° (Figure 5). In contrast, for a larger value of α, the central angular momentum of the simulation with θ = 0 is greater than that with θ = 90° (Figure 6). These results are consistent with previous studies that employed a large value of α (Matsumoto & Tomisaka 2004) and a small α (Hennebelle & Ciardi 2009; Joos et al. 2012). The results become simpler once magnetic diffusion has been incorporated into the simulation. The central angular momentum decreases as θ increases, regardless of the value of α (Figure 6).

3.3. Angular Momentum Evolution of the Fluid Element in the Central Region

This subsection investigates the underlying mechanism that causes the diversity of the evolution of the angular momentum described above. In the previous subsection and in most earlier studies, the temporal evolution of the central angular momentum of a collapsing cloud core was investigated by specifying threshold densities for the "central region" (e.g., Joos et al. 2012; Tsukamoto et al. 2015b) or a criterion that determines the "disk" (Masson et al. 2016). However, these kinds of analyses make it difficult to identify the mechanism causing the diversity of the evolution of the angular momentum, because the resultant angular momentum of the central region is determined by the supply of angular momentum provided by the envelope accretion (the angular momentum of the accretion flow may have changed from its initial value already) and by the removal of the angular momentum from the central region by the magnetic field. Thus, it is unclear when and how the angular momentum of the gas changes. We therefore investigate the evolution of the angular momentum of each fluid element during cloud core collapse.

The following three subsections show that three different mechanisms affect the evolution of the angular momentum evolution. The first is anisotropic accretion caused by radial magnetic tension. The second is magnetic braking during the isothermal collapse phase. The third is magnetic braking in the first-core or newborn disk, which plays a significant role in the ideal MHD simulations.

3.3.1. Initial Locations of the Gas Falling into the Central Region

Figure 8 shows the initial distribution of the fluid elements that have fallen into the central region by the end of the ideal MHD simulations (by the epoch when the central density has reached ${\rho }_{{\rm{c}}}\sim {10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ). Following the criterion used in the previous section, we define the central region as the region where $\rho \gt {10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ . At the end of each simulation run, 10⁴ SPH particles were sampled from the central region. We then traced their trajectories back and determined their initial positions.

**Figure 8.** Initial distribution of the SPH particles that have fallen into the central region (with density greater than ${10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ at the end of the simulation), obtained from the ideal MHD simulations (models I1–I9). The color scale shows the specific angular momentum of each particle in the initial state. The initial angular momentum of the core is parallel to the z-axis, and the amount of specific angular momentum of each particle is proportional to (x² + y²).
Download figure:
Standard image High-resolution image

The figure shows that the distributions of the initial fluid element locations are not spherically symmetric but prolate, with the major axis running nearly parallel to the initial magnetic field. During cloud core collapse, the gas tends to move along the field line by the Lorentz force, thereby giving the distributions of the prolate shape.

The color of each fluid element represents the initial angular momentum. Because we assumed that the rotation axis of the initial cloud core was parallel to the z-axis, the angular momentum of the fluid elements is proportional to (x² + y²), i.e., the square of the distance from the z-axis. Because the major axis of the prolate distribution is found to be parallel to the initial magnetic field, the fluid elements that initially had a small angular momentum selectively accreted onto the central region when θ = 0°, whereas the fluid elements that initially have a large angular momentum selectively accreted onto the central region when θ = 90°. If we only investigate the central angular momentum, this nonspherical symmetric accretion apparently causes strong and weak magnetic braking in the cases with, respectively, θ = 0° and 90°. Note that this effect is "apparent" in the sense that the fluid elements with a large (for θ = 0°) or small (for θ = 90°) angular momentum eventually accrete to the central region in the subsequent evolution, though their timing is delayed.

By comparing the panels with the same θ values (e.g., panels (a), (b), and (c)), we find that the ratio of the major axis of the particle distribution decreases as α decreases and that the initial distribution becomes more spherical. A core with a small α corresponds to a more unstable cloud core. A stronger gravitational force is exerted on each fluid element, and the Lorentz force becomes relatively weak. As a result, the initial distribution is more spherical in simulations with a small α.

3.3.2. Evolution of the Angular Momentum of the Fluid Elements in the Central Region

Figure 9 shows the temporal evolution of the mean specific angular momenta of the fluid elements that have fallen into the central region at the end of the ideal MHD simulations (the particles shown in Figure 8). The horizontal axis represents the median density of the fluid elements, which describes the time evolution. We use the median density of the sampled particles for the following reasons. Using Figures 9 and 10, we want to determine the typical density where the most sampled particles lose their angular momentum. However, we found that the central density does not correctly represent the typical density. We tried to find a good quantity for our purpose and found the median density to be suitable for indicating the density where most particles lose their angular momentum. The reason for the initial increase in angular momentum as θ increases (indicated with arrows (i) in Figure 8)) is the nonspherical accretion discussed in the previous subsection.

**Figure 10.** Temporal evolution of the mean specific angular momentum of the fluid elements that have fallen into the central region (with density greater than ${10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ at the end of the non-ideal MHD simulations). The horizontal axis represents the median density of the fluid elements, which describes the temporal evolution. The vertical axis shows the mean specific angular momentum of the fluid elements. The arrows indicate the positions where (i) the initial angular momentum, (ii) magnetic braking in the isothermal collapse phase, and (iii) magnetic braking in the first core and disk are found.
Download figure:
Standard image High-resolution image

Figure 9 indicates that the mean specific angular momentum decreases in the two different density ranges. A decrease occurs in $\rho \lesssim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , i.e., in the isothermal collapse phase (indicated with arrows (ii)). At this stage, the angular momentum decreases more rapidly in the models with a larger θ. This suggests that the magnetic braking in the isothermal phase is strong for a greater θ. We will come back to this point in the next subsection.

The other decrease occurs in $\rho \gtrsim {10}^{-11}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , i.e., in the first core or disk around it (indicated with arrows (iii)). In this phase, the gas is supported by the pressure or rotation (when the disk around the first core is formed; see Tomida et al. 2015; Tsukamoto et al. 2015b). The two-step decrease of the angular momentum was observed in the pioneering 2D ideal MHD simulation by Tomisaka (2000). He showed that approximately 70% of the angular momentum is extracted during the isothermal collapse phase and that the remaining angular momentum is extracted during the adiabatic collapse phase. The results of our ideal MHD simulations with θ = 0° are consistent with the results obtained in Tomisaka (2000).

Figure 9 also shows that the decreases in angular momentum in both the isothermal and adiabatic phases are less significant in a cloud core with a small α. Because gravity is strong in the core with a small α and freefall time is short, there is not enough time for the magnetic field to extract the angular momentum from the fluid elements. The results again suggest that α (or the degree of gravitational instability of the initial core) is an important parameter for investigating the effect of magnetic braking in collapsing cloud cores.

Figure 10 shows the evolution of the mean specific angular momentum for the central fluid elements in the non-ideal MHD simulations. As in the case of Figure 9, 10⁴ SPH particles were sampled from the region with $\rho \gt {10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ at the end of the simulation, and their trajectories were traced back to determine their initial positions. Note that the sampled particles are not the same in the ideal and non-ideal MHD simulations for given values of α and θ, despite starting from the same initial cloud core. This is because the structures of the central region at the end of the simulations are different as a result of inclusion of the magnetic diffusion.

Figure 10 shows that the decrease in angular momentum in the isothermal collapse phase ( $\rho \lesssim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) also occurs in non-ideal MHD simulations (indicated with arrows (ii)) and that the mean angular momentum in the θ = 90° model becomes smaller than that for θ = 0° in the isothermal phase in all simulations. The evolution of the angular momentum in the isothermal phase is similar to that obtained with the ideal MHD approximation because magnetic diffusion is effective only for $\rho \gtrsim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ .

Unlike the results of the ideal MHD simulations, the secondary decrease in the angular momentum (indicated with arrows (iii)) is weak. The ohmic and ambipolar diffusion efficiently removes the magnetic flux from the region with $\rho \gtrsim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , and the magnetic field is weaker in the central region (see Figures 3 and 4 of Tsukamoto et al. 2015b) than in the ideal MHD simulations. Furthermore, the gas and the magnetic field are almost decoupled in the central region and the toroidal magnetic field generated by the gas rotation, which causes the magnetic torque to also be weak. As a result, the efficiency of magnetic braking is reduced in the high-density region compared to the ideal MHD simulations.

3.3.3. Evolution of the Angular Momentum of the Fluid Elements Initially in the Spherical Region

The initial distributions of the fluid elements discussed in the previous section differ between models because the fluid elements are sampled at the end of the simulations (Figure 8). Thus, although we have shown that the angular momentum significantly decreases in the isothermal collapse phase in the simulations with θ = 90°, it remains debatable whether the magnetic braking in the isothermal collapse phase is stronger in the perpendicular than in the parallel configuration. This is because the initial distribution discussed in the previous section is elongated along the initial magnetic field and has a longer lever arm in the core with a greater θ. The long and small lever arm may enhance magnetic braking when θ = 90° and suppress it when θ = 0° in the previous section.

To clarify the issue decisively, we investigated the evolution of the angular momentum of the fluid elements of a sphere that is selected from the initial cloud core. Let the sphere be parameterized by the enclosed mass,

$\begin{eqnarray}&&{M}_{{\rm{e}}}(r)={\int }_{0}^{r}\rho (r)d{\boldsymbol{V}}.\end{eqnarray} \tag{ 7 }$

We chose a sphere with ${M}_{{\rm{e}}}(r)=5\times {10}^{-2}\,{M}_{\odot }$ , i.e., the innermost region of the entire cloud core. Then 10⁴ SPH particles were sampled from this sphere at the beginning of the simulation, and their trajectories were traced. By this sampling procedure, the fluid elements of the simulations with the same α are the same, and we can investigate the difference of the magnetic braking efficiency purely caused by misalignment.

Figure 11 shows the evolution of the mean angular momentum of the sphere in the isothermal collapse phase. The results of both the ideal (lines) and non-ideal (symbols) MHD simulations are shown. The figure clearly indicates that the angular momentum becomes smaller in the simulations with θ = 90° (dotted lines) than with θ = 0° (solid lines) during the isothermal collapse phase. The initial distributions of fluid elements in each panel are identical; hence, the difference in the angular momentum is caused purely by the difference of the magnetic braking efficiency that originates in misalignment. Therefore, we conclude that magnetic braking in the isothermal collapse phase is more effective in the perpendicular than in the parallel configuration.

**Figure 11.** Temporal evolution of the mean specific angular momentum of the fluid elements in a sphere with ${M}_{e}(r)=5\times {10}^{-2}\,{M}_{\odot }$ . The results of the ideal (lines) and non-ideal (symbols) MHD simulations are shown (models I1–I9 and NI1–NI9, respectively). The dot-dashed lines show the results of the hydrodynamics simulations. The horizontal axis shows the median density of the fluid elements. The vertical axis shows the mean specific angular momentum of the fluid elements.
Download figure:
Standard image High-resolution image

Another important finding is that, with smaller α values, magnetic braking in the isothermal collapse phase becomes weak, and the difference in the angular momenta between simulations with different θ values is small. Therefore, in the simulations with small α values, magnetic braking in the high-density region plays a crucial role in determining whether the central angular momentum of θ = 0° is greater than that of θ = 90°. The figure also shows that the evolution of the angular momentum is almost identical in the isothermal collapse phase between the ideal and non-ideal MHD simulations.

The temporal evolution of the angular momentum of the hydrodynamical simulations (more precisely, ideal MHD simulations with a very weak magnetic field of μ = 10⁴) is plotted with dot-dashed lines. The lines show that the angular momentum is almost constant during the isothermal collapse phase. This confirms that the decrease of the angular momentum observed in MHD simulations is indeed caused by the magnetic field.

4. Discussion

4.1. Three Mechanisms for the Evolution of the Angular Momentum during the Cloud Core Collapse and Their Characteristics

In Section 3.3, we showed that three mechanisms affect the angular momentum evolution of the central region, i.e.,

(i)
the selective accretion of the rapidly (for the core with θ = 90°) or slowly (for the core with θ = 0°) rotating fluid elements to the central region;
(ii)
magnetic braking in the isothermal collapse phase; and
(iii)
magnetic braking in the first core (adiabatic collapse phase) or the disk.

We summarize the characteristics of these mechanisms and their dependence on the initial conditions and magnetic diffusions.

We have shown that magnetic braking in the isothermal collapse phase is stronger in the cores with θ = 90° than with θ = 0° by following the evolution of the angular momentum of the fluid elements. This result is consistent with the analytic argument by Mouschovias (1985) and the simulation result by Matsumoto & Tomisaka (2004). The enhancement of magnetic braking by flared magnetic field geometry that was proposed as a cause for stronger magnetic braking in the case of θ = 0° in Joos et al. (2012) does not seem to be effective in the isothermal collapse phase.

Joos et al. (2012) pointed out that the timescale of magnetic braking of a core with θ = 0° and a flared magnetic field geometry becomes smaller than that with θ = 90°. According to their estimate, the magnetic braking timescale of the core with θ = 0° and a field geometry of panel (c) of Figure 12 is smaller than that with θ = 90° and a field geometry of panel (b) of Figure 12 by a factor of (R_c/R_core) (≪1) times, where R_c and R_core are the radius of the central region and of the initial core, respectively. This argument seems to contradict our conclusions that magnetic braking during isothermal collapse phase tends to be strong in the cores with θ = 90°. An explanation is therefore required. In the Appendix, we compare the timescales for magnetic braking with various field geometries and show that, even with a flared magnetic field geometry, the timescale for magnetic braking for the core with θ = 0° is not always smaller than that with θ = 90°. Thus, our results do not contradict the analytical estimate. Rather, the analytical estimates have large uncertainties, and even the qualitative result depends on the assumptions made.

The efficiency of magnetic braking in the isothermal collapse phase decreases as α decreases (Figures 9 and 11). In this study, we parameterized the initial cloud core with $\alpha \equiv {E}_{\mathrm{therm}}/{E}_{\mathrm{grav}}$ . Because the initial temperature is fixed to 10 K, a small α corresponds to a large gravitational energy and short freefall timescale. The magnetic field does not then have enough time to change the angular momentum of the fluid elements during the isothermal collapse phase, which in turn makes the magnetic braking become weak in the core with small α.

The selective accretion apparently weakens magnetic braking in the simulations with θ = 90° and strengthens it with θ = 0°. The radial magnetic tension suppresses mass accretion from the direction perpendicular to the field line and selectively enhances mass accretion parallel to the field line in the early evolution phase. As a result, the fluid elements with a greater initial angular momentum selectively accrete onto the central region in the simulations with θ = 90°, while those with a smaller angular momentum selectively accrete with θ = 0°. The effect of selective accretion is less significant in a core with a smaller α (Figure 8).

Selective accretion may serve as an additional mechanism for increasing the central angular momentum in a core with θ = 90° in its early evolution phase. However, this effect is "apparent" in the sense that fluid elements with a smaller angular momentum eventually accrete to the central region in the subsequent evolution. In contrast, the disk is expected to grow in the later evolutionally phase in the core with θ = 0° because the gas with a larger angular momentum eventually accretes to the center. Indeed, Machida et al. (2011) reported disk growth in the later evolutionary phase with long-term simulations up to 10⁵ yr after the protostar formation. They suggested that the growth of the disk in the late phase is caused by the decrease in the magnetic braking efficiency brought about by the depletion of the envelope above the circumstellar disk. However, we suggest that the delay in the accretion of fluid elements with a large angular momentum serves as another mechanism for the rapid increase of the disk size in the late phase.

Of the three mechanisms, magnetic braking in the high-density region (in the first core and disk) may depend weakly on the gravitational energy of the initial cloud core because the gas in the high-density region is supported by the pressure gradient force or the centrifugal force, and because the freefall time of the initial cloud core is not the characteristic timescale. Owing to the pressure and rotation support, the magnetic field has much more time than the local freefall time to extract the gas angular momentum. Because the magnetic field fans out around the central region, the magnetic braking in the high-density region of the cores with θ = 0° may be enhanced by the hourglass-like magnetic field geometry as indicated by Equations (23) and (25). This, along with nonspherical accretion, may cause the flipping of the central angular momentum between the core with θ = 0° and the core with θ = 90° in ideal MHD simulations of α = 0.2.

4.2. Importance of the Non-ideal Effect

The efficiency of magnetic braking in the high-density region depends strongly on whether magnetic diffusions are taken into account. In the ideal MHD simulations, the significant decrease in the angular momentum is observed in the high-density region (Figure 9). Once magnetic diffusions have been taken into account, however, the decrease in the angular momentum becomes moderate (see arrows (iii) in Figure 10). The ambipolar and ohmic diffusions are effective in the region with $\rho \gtrsim {10}^{-13}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and in $\rho \gtrsim {10}^{-11}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ , respectively (Dapp et al. 2012; Tomida et al. 2015; Tsukamoto et al. 2015b; Masson et al. 2016; Wurster et al. 2016), and the magnetic flux is removed from the central region. Furthermore, induction of the toroidal magnetic field by gas rotation is also suppressed. As a result, the toroidal magnetic tension in the high-density region is weak in the non-ideal MHD simulations. The simulations of disk formation with the ideal MHD simulations may exaggerate the impact of magnetic braking in the high-density region. Our results show that the difference between the ideal and non-ideal MHD simulations is significant and suggest that the inclusion of magnetic diffusion is crucial for investigating circumstellar disk formation and evolution.

The effect of the magnetic diffusions depends sensitively on the resistivity models (e.g., cosmic-ray ionization rate) and the dust models (Nishi et al. 1991; Dapp et al. 2012; Dzyurkevich et al. 2017; Zhao et al. 2018a). Their difference must therefore be borne in mind when comparing the previous studies with non-ideal effects (e.g., Tomida et al. 2015; Tsukamoto et al. 2015b; Masson et al. 2016; Wurster et al. 2016, 2018a; Zhao et al. 2016, 2018b; Vaytet et al. 2018) quantitatively, although our results agree qualitatively with those earlier studies. A thorough systematic study of the effect of the variety of the resistivity model remains to be done (see, however, the impact of the difference in cosmic-ray ionization rate in Wurster et al. 2018a, 2018b). Furthermore, it would be worthwhile to investigate the impact of dust growth in the disk, which may happen even in the earliest phase of the disk evolution (Tsukamoto et al. 2017b). We believe that the detailed study of the resistivity models is an important subject for future work.

The present study neglects the Hall effect. Its impact on the evolution of angular momentum in misaligned cloud cores is studied in detail in Tsukamoto et al. (2017a), which considers a uniform cloud core with α = 0.3 and a uniform magnetic field with μ = 4. The Hall effect induces rotation in the left-handed screw direction of magnetic field during the isothermal collapse phase. As a result, the angular momentum of the central region (e.g., $\rho \gt {10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ ) of the core with θ = 180° is an order of magnitude greater that that with θ = 0° in Tsukamoto et al. (2017a). The angular momenta of these two cases become identical when the Hall effect is neglected. The angular momentum of the core with θ = 90° has an intermediate value between these two cases. Tsukamoto et al. (2017a), however, only investigated a single value of α, and it is therefore unclear how the impact of Hall effect changes in the core with different values of α. Magnetic braking in the isothermal collapse phase can be expected to become more significant with a greater α. Therefore, the quantitative difference between cores with different θ values may depend on α. We intend to investigate this issue in future works.

4.3. Comparison with Previous Studies

Matsumoto & Tomisaka (2004) reported that magnetic braking is more efficient in a core with θ = 90° than with θ = 0° (Figure 12 of Matsumoto & Tomisaka 2004). Their simulations adopted the ideal MHD approximation and modeled the cloud core as a Bonnor–Ebert sphere with a uniform magnetic field. The parameter α of the core was α = 0.5. The simulations were terminated at the protostar formation epoch. Although the density profile is different from our initial condition, their result is consistent with ours, obtained from ideal MHD simulations with α = 0.6 and α = 0.4. With such large α values, magnetic braking in the isothermal collapse phase plays a crucial role in the evolution of angular momentum.

Hennebelle & Ciardi (2009) and Joos et al. (2012) reported that magnetic braking is more efficient in a core with θ = 0° than with θ = 90°, in apparent contradiction with results by Matsumoto & Tomisaka (2004). Note, however, that the difference in the results of Hennebelle & Ciardi (2009), Joos et al. (2012), and Matsumoto & Tomisaka (2004) is not contradictory but originates from different samplings of the parameter space, as we discuss below. Hennebelle & Ciardi (2009) and Joos et al. (2012) also adopted the ideal MHD approximation, and the cloud core was modeled as a gas sphere with a Plummer-like density profile. The magnetic field strength in Hennebelle & Ciardi (2009) was set to be proportional to the column density, whereas no mention of the configuration of the magnetic field was given in Joos et al. (2012). Their cores have α = 0.25. Using a barotropic EOS in which the second collapse is artificially suppressed, Joos et al. (2012) followed the simulation for several thousand years after protostar formation. Figure 4 in Joos et al. (2012) shows that the angular momentum of the central region of the core with θ = 0° is smaller than that with θ = 90°. Although the density profile is different from our initial condition, this result is consistent with ours, obtained with ideal MHD simulations with α = 0.2. With this small α value, magnetic braking in the isothermal collapse phase is weak, and the evolution of the central angular momenta depends more strongly on the magnetic braking efficiency in the high-density region.

Indeed, the results by Joos et al. (2012) support our conclusion that the evolution of the central angular momentum in the core with a small α value is determined by magnetic braking in the high-density region. Figure 4 in Joos et al. (2012) shows that the difference in the mean angular momenta between models with different θ values is only apparent in the high-density region ${n}_{{\rm{c}}}\gt {10}^{10}{\mathrm{cm}}^{-3}$ and that the mean angular momentum, in the range including the low-density region ( ${n}_{{\rm{c}}}\gt {10}^{8}{\mathrm{cm}}^{-3}$ ), is almost independent of θ. This result suggests that magnetic braking in the isothermal collapse phase does not introduce a significant difference between the simulations with different θ values and that the difference is caused by magnetic braking in the high-density region.

We therefore conclude that the apparent discrepancy between Matsumoto & Tomisaka (2004) and Joos et al. (2012) is due to the magnetic braking efficiency in the isothermal collapse phase. Matsumoto & Tomisaka (2004) employed a relatively large α value, and a significant fraction of the angular momentum is removed during the isothermal collapse phase in their simulations with θ = 90°. Whether the central angular momenta of θ = 0° is greater than that of θ = 90° is determined in this phase. This is inferred also from Figure 12 of Matsumoto & Tomisaka (2004), which shows that the difference in the central angular momenta is introduced in the isothermal collapse phase. In contrast, Hennebelle & Ciardi (2009) and Joos et al. (2012) adopted a smaller α value than Matsumoto & Tomisaka (2004). A smaller α value makes the magnetic braking in the isothermal collapse phase less significant, and the difference in the angular momenta at the end of the isothermal collapse phase between models with θ = 0° and θ = 90° becomes small (see Figures 9 and 11). Because magnetic braking in the high-density region is significant in the ideal MHD simulation and is possibly stronger in the core with θ = 0° than with θ = 90° owing to the hourglass-like magnetic field (see the Appendix), the central angular momentum in the core with θ = 0° is smaller than that with θ = 90°. The selective accretion may serve as an additional mechanism to make the angular momentum of θ = 0 become smaller than that of θ = 90.

We note, however, that the magnetic braking in the high-density region is strongly affected by whether the magnetic diffusions are included. Figures 5, 6, 9, and 10 show significant differences in the evolution of the angular momentum in the high-density region.

Masson et al. (2016) reported that the small angular misalignment of ${\boldsymbol{B}}$ and ${{\boldsymbol{J}}}_{\mathrm{ang}}$ (θ < 40°) does not change the disk angular momentum noticeably. They considered ohmic and ambipolar diffusions and modeled the cloud core as a uniform gas sphere with a constant magnetic field and α = 0.25. They adopted a barotropic EOS and conducted simulations for several thousand years after the formation of the first core. Figure 13 in Masson et al. (2016) shows that the evolution of angular momentum of the disk is almost identical for θ = 0° and θ = 40°. Their results are consistent with ours, shown in Figure 6, because there is only a small difference between the cases with θ = 0° and 45°.

Once magnetic diffusion is taken into account, magnetic braking in the high-density region is suppressed (arrows (iii), in Figure 10) and the effect of magnetic braking in the isothermal phase is preserved in the central region. As a result, the central angular momentum decreases as θ increases, independently of α (Figure 6). We conclude that this is the primary reason why a significant difference in angular momentum is not observed in Masson et al. (2016), where α = 0.25 and θ < 40°.

In consideration of previous studies, our results are inconsistent with those of Li et al. (2013), who considered a core with uniform density with α = 0.7 and employed ideal MHD simulations. They reported that circumstellar disks form more easily in cores with θ = 90° than with θ = 0° even with such a large value for α. This contradicts our simulation results and the results of Matsumoto & Tomisaka (2004).

The reason for this discrepancy is unclear. One possible reason is the different treatment of the inner boundary. Li et al. (2013) used an open inner boundary and removed the gas from the system. The gas angular momentum was also removed simultaneously. This treatment is different from that used in other studies (Matsumoto & Tomisaka 2004; Hennebelle & Ciardi 2009; Joos et al. 2012; Masson et al. 2016), in which an open inner boundary was not used and the gas accumulated in the central region. The first core and the disk around it serve as a reservoir for the angular momentum. Thus, the angular momentum also accumulates around the center. In contrast, the use of an open inner boundary suppresses such an accumulation of angular momentum.

Although we have reproduced almost all the qualitative results reported in previous studies, some inconsistencies remain, as discussed above. These may result from differences in other simulation configurations, e.g., boundary conditions, initial cloud cores, the EOS, or numerical schemes. We think that the analysis of the evolution of angular momentum of the fluid elements provides precise information on magnetic braking. For future studies of this subject, we suggest analyzing the evolution of the angular momentum of Lagrangian fluid elements as described in this paper.

5. Conclusions

Our non-ideal MHD simulations have revealed that circumstellar disks are formed more easily in a core where the rotation axis is aligned with the magnetic field direction, independently of the gravitational stability parameter α. This is because magnetic braking in the isothermal collapse phase is stronger in a core with θ = 90° than with θ = 0°. Also, magnetic diffusion suppresses magnetic braking in the first core and circumstellar disk. The role of magnetic diffusions has shown to be crucial, and ideal MHD simulations may exaggerate magnetic braking in circumstellar disks.

We thank Tomoyuki Hanawa and Tomoaki Matsumoto for fruitful discussion. We also thank an anonymous referee for helpful comments that greatly improved this paper. We thank K. Tomida and Y. Hori for providing us with their equation of state, which was used in Tomida et al. (2013). The computations were performed on a parallel computer, XC40/XC50 system at the Center for Computational Astrophysics of National Astronomical Observatory of Japan. This work is supported by JSPS KAKENHI grant nos. 17KK0096, 17K05387, 17H06360, 17H02869, and 18K13581.

Appendix: Analytical Discussion of the Magnetic Braking

In this appendix, we derive the magnetic braking timescales based on the analytical model of Mouschovias & Paleologou (1979, 1980) and Mouschovias (1985). We particularly emphasize that a flared magnetic field geometry does not always cause the order reversal of the magnetic braking timescale for the parallel and perpendicular configurations. Mouschovias & Paleologou (1979, 1980) showed that the magnetic braking timescale t_b of a central region with a moment of inertia I_c can be estimated as the timescale over which Alfvén waves sweep through an amount of gas in the outer envelope with a moment of inertia I_ext(t_b) equal to I_c. This condition is given by

$\begin{eqnarray}&&{I}_{\mathrm{ext}}({t}_{{\rm{b}}})={I}_{{\rm{c}}}.\end{eqnarray} \tag{ 8 }$

This then determines the magnetic braking timescale if the structures of the central region, outer envelope, and magnetic field are specified.

The magnetic braking timescale for a simple parallel configuration ( ${{\boldsymbol{J}}}_{\mathrm{ang}}\parallel {\boldsymbol{B}}$ or θ = 0°) can be calculated as follows. The central collapsing region is modeled as a uniform cylinder with a density ${\rho }_{\mathrm{cyl}}$ , radius R_cyl, and height $2{H}_{\mathrm{cyl}}$ , where a uniform magnetic field runs parallel to the rotation axis and the density of the outer envelope ρ_ext is assumed to be constant (panel (a) of Figure 12). The moments of inertia of the central cylinder and outer envelope are given by ${I}_{{\rm{c}}}=\pi {\rho }_{\mathrm{cyl}}{R}_{\mathrm{cyl}}^{4}{H}_{\mathrm{cyl}}$ and ${I}_{\mathrm{ext}}({t}_{{\rm{b}}})=\pi {\rho }_{\mathrm{ext}}{R}_{\mathrm{cyl}}^{4}{v}_{{\rm{A}}}{t}_{{\rm{b}}}$ , respectively, where v_A denotes the Alfvén velocity in the outer envelope and v_At_b corresponds to the height of the swept outer envelope. By solving Equation (8), t_b is calculated to be

$\begin{eqnarray}&&{t}_{{\rm{b}},\parallel }=\displaystyle \frac{{\rho }_{\mathrm{cyl}}}{{\rho }_{\mathrm{ext}}}\displaystyle \frac{{H}_{\mathrm{cyl}}}{{v}_{{\rm{A}}}}\end{eqnarray} \tag{ 9 }$

(Mouschovias 1985). Using the mass $M=2\pi {\rho }_{\mathrm{cyl}}{R}_{\mathrm{cyl}}^{2}{H}_{\mathrm{cyl}}$ and the magnetic flux ${\rm{\Phi }}=\pi {R}_{\mathrm{cyl}}^{2}B$ of the cylinder, we can reduce t_b,∥ to

$\begin{eqnarray}&&{t}_{{\rm{b}},\parallel }={\left(\displaystyle \frac{\pi }{{\rho }_{\mathrm{ext}}}\right)}^{1/2}\displaystyle \frac{M}{{\rm{\Phi }}}.\end{eqnarray} \tag{ 10 }$

This shows that the magnetic braking timescale in the simple geometry depends only on the mass-to-flux ratio of the central region and the density of the outer envelope.

**Figure 12.** Schematic diagrams of the magnetic field geometries discussed in the Appendix. Panels (a), (b), (c), and (d) show the geometry for Equations (9), (11), (15), and (23), respectively.
Download figure:
Standard image High-resolution image

The magnetic braking timescale in a simple geometry with a perpendicular configuration ( ${{\boldsymbol{J}}}_{\mathrm{ang}}\perp {\boldsymbol{B}}$ or θ = 90°) is calculated as follows. The central collapsing region is modeled in the same way as with the parallel configuration, using the parameters ${\rho }_{\mathrm{cyl}}$ , R_cyl, and $2{H}_{\mathrm{cyl}}$ . The magnetic field is assumed to be perpendicular to the rotation axis and axisymmetrical (panel (b) of Figure 12) where $B(r)\propto {r}^{-1}$ . The density of the outer envelope ρ_ext is again assumed to be constant (panel (b) of Figure 12).

The moments of inertia of the central cylinder and outer envelope are given by ${I}_{{\rm{c}}}=4\pi {\rho }_{\mathrm{cyl}}{H}_{\mathrm{cyl}}{R}_{\mathrm{cyl}}^{4}$ and ${I}_{\mathrm{ext}}\,=4\pi {\rho }_{\mathrm{ext}}{H}_{\mathrm{cyl}}({R}_{\mathrm{ext}}^{4}-{R}_{\mathrm{cyl}}^{4})$ . Using the condition ${I}_{{\rm{c}}}={I}_{\mathrm{ext}}$ , the magnetic braking timescale for the perpendicular configuration is calculated to be

$\begin{eqnarray}\begin{array}{rcl}{t}_{{\rm{b}},\perp } & = & {\int }_{{R}_{\mathrm{cyl}}}^{{R}_{\mathrm{ext}}}\displaystyle \frac{{dr}}{{v}_{{\rm{A}}}}=\displaystyle \frac{1}{2}\left({\left(\displaystyle \frac{{\rho }_{\mathrm{cyl}}}{{\rho }_{\mathrm{ext}}}+1\right)}^{1/2}-1\right)\displaystyle \frac{{R}_{\mathrm{cyl}}}{{v}_{{\rm{A}}}({R}_{\mathrm{cyl}})}\\ & & \sim \,\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{\rho }_{\mathrm{cyl}}}{{\rho }_{\mathrm{ext}}}\right)}^{1/2}\displaystyle \frac{{R}_{\mathrm{cyl}}}{{v}_{{\rm{A}}}({R}_{\mathrm{cyl}})}\end{array}\end{eqnarray} \tag{ 11 }$

(Mouschovias 1985), where we assumed ${\rho }_{\mathrm{cyl}}/{\rho }_{\mathrm{ext}}\gg 1$ . In addition, t_b,⊥ is reduced to

$\begin{eqnarray}&&{t}_{{\rm{b}},\perp }\sim 2{\left(\displaystyle \frac{\pi }{{\rho }_{\mathrm{cyl}}}\right)}^{\tfrac{1}{2}}\displaystyle \frac{M}{{\rm{\Phi }}}.\end{eqnarray} \tag{ 12 }$

Note that this formula is likely to overestimate the timescale (or to underestimate the magnetic braking efficiency) because, in the magnetic field configuration, the Alfvén velocity decreases as ${v}_{{\rm{A}}}\propto {r}^{-1}$ and a long lever arm cannot be obtained. If the magnetic field has the anisotropic configuration, the Alfvén wave can propagate further and the longer lever arm decreases the braking timescale.

The ratio of the magnetic braking timescale for the perpendicular and the parallel configurations is

$\begin{eqnarray}&&\displaystyle \frac{{t}_{{\rm{b}},\perp }}{{t}_{{\rm{b}},\parallel }}\sim {\left(\displaystyle \frac{{\rho }_{\mathrm{ext}}}{{\rho }_{\mathrm{cyl}}}\right)}^{1/2}.\end{eqnarray} \tag{ 13 }$

Thus, the timescale for perpendicular configuration is shorter than that for the parallel one because ${\rho }_{\mathrm{cyl}}\gg {\rho }_{\mathrm{ext}}$ and the magnetic braking in the former is stronger than in the latter. This is the straightforward conclusion derived for the simplest geometry cases.

However, the simple geometry of the parallel configuration may be inappropriate for a model of a collapsing cloud core. The magnetic field is suggested to have an hourglass-like shape as a result of radial dragging according to theoretical work (e.g., Tomisaka 2002; Allen et al. 2003; Kunz & Mouschovias 2010; Tsukamoto et al. 2015b; Vaytet et al. 2018) and as has been confirmed observationally in YSOs (Girart et al. 2006). In the hourglass configuration, the magnetic field fans out vertically, which is inconsistent with the simple geometry of the parallel configuration.

Mouschovias (1985) derived the magnetic braking timescale in which the magnetic flux tube expands quickly in an infinitely thin transition layer to a radius R_ext just above the central cylinder (panel (c) of Figure 12). If we ignore the moment of inertia of the transition layer, I_ext(t_b) is given by

$\begin{eqnarray}&&{I}_{\mathrm{ext}}({t}_{{\rm{b}}})=\pi {\rho }_{\mathrm{ext}}{R}_{\mathrm{ext}}^{4}{v}_{{\rm{A}}}{t}_{{\rm{b}}}.\end{eqnarray} \tag{ 14 }$

Using ${I}_{{\rm{c}}}=\pi {\rho }_{\mathrm{cyl}}{R}_{\mathrm{cyl}}^{4}{H}_{\mathrm{cyl}}$ and Equation (14), we obtain the magnetic braking timescale of the disk with an hourglass magnetic field geometry (Mouschovias 1985) as

$\begin{eqnarray}&&{t}_{{\rm{b}},{\rm{f}}}={\left(\displaystyle \frac{\pi }{{\rho }_{\mathrm{ext}}}\right)}^{1/2}\left(\displaystyle \frac{M}{{\rm{\Phi }}}\right){\left(\displaystyle \frac{{R}_{\mathrm{cyl}}}{{R}_{\mathrm{ext}}}\right)}^{2},\end{eqnarray} \tag{ 15 }$

where

$\begin{eqnarray}&&{B}_{\mathrm{ext}}={\left(\displaystyle \frac{{R}_{\mathrm{cyl}}}{{R}_{\mathrm{ext}}}\right)}^{2}{B}_{\mathrm{cyl}}\end{eqnarray} \tag{ 16 }$

is assumed because of the conservation of magnetic flux. According to Equation (15), the magnetic braking timescale can be shorter than t_b,∥ in Equation (10) because of the factor (R_cyl/R_ext)² (<1).

Joos et al. (2012) argued that the magnetic braking timescale becomes smaller in the parallel than perpendicular configuration by virtue of Equation (15). They assumed R_ext = R_core, where R_core is the initial core radius, and ρ(r) ∝ r⁻². Then, assuming that the mean external density and the central density can be approximated in terms of the volume-averaged means as, respectively, ${\rho }_{\mathrm{cyl}}\sim \rho ({R}_{\mathrm{cyl}})$ and ${\rho }_{\mathrm{ext}}\sim \rho ({R}_{\mathrm{ext}})$ , the ratio of the timescales is estimated to be

$\begin{eqnarray}&&\displaystyle \frac{{t}_{{\rm{b}},\perp }}{{t}_{{\rm{b}},{\rm{f}}}}\sim \displaystyle \frac{{R}_{\mathrm{core}}}{{R}_{\mathrm{cyl}}}.\end{eqnarray} \tag{ 17 }$

Using this relation, they concluded that the braking timescale in flared parallel configuration is shorter than in the perpendicular configuration because R_cyl < R_core.

However, it is uncertain whether the assumption of an infinitely thin magnetic flux tube expansion layer (transition layer) is valid in realistic situations, given that such a layer causes a strong dependence of the braking timescale on the ratio of the radii t_b,f ∝ (R_cyl/R_ext)². This assumption corresponds to a discontinuous decrease of poloidal magnetic field strength (Equation (16)) at the vertical height z = H_cyl. The poloidal magnetic field in a collapsing cloud core, however, may change in a power-law fashion. Furthermore, the analytic estimate ignores the propagation time of the Alfvén wave in the expansion layer, which would not be negligible when R_cyl ≪ R_ext.

We should also note that the validity of the assumption of R_ext = R_core is unclear, because R_ext is the radius of the external flux tube and is not related to that of the core R_core, and similarly because the radius in the cylindrical coordinates (or distance from the vertical axis) is not related to the radius in the spherical coordinates (or distance from the center). R_ext can be R_ext ≪ R_core even at the vertical height of z = R_core unless all magnetic flux of the core is concentrated in the central region, which is unlikely in the early phase of disk formation. Therefore, Equation (15) should be regarded as the lower limit of the magnetic braking timescale in the parallel configuration.

To demonstrate clearly that the discontinuous expansion of magnetic flux tubes and the assumptions described above lead to the flip of the magnetic braking timescale, we calculate the magnetic braking timescale with magnetic flux tubes that expand in a power-law fashion. Our result will show that whether the timescale becomes shorter or not depends on the power-law index and external density in this case.

We consider the configuration where the radius of the magnetic flux tube obeys

$\begin{eqnarray}&&{R}_{\mathrm{ext}}(z)={R}_{\mathrm{cyl}}{\left(\displaystyle \frac{z}{H}\right)}^{\beta }\end{eqnarray} \tag{ 18 }$

(see panel (d) of Figure 12). We consider the case of β > 0 only, i.e., where the flux tube is expanding. The magnetic field strength is approximated to

$\begin{eqnarray}&&B(z)\sim {B}_{\mathrm{cyl}}\displaystyle \frac{{R}_{\mathrm{cyl}}^{2}}{{R}_{\mathrm{ext}}{(z)}^{2}}={B}_{\mathrm{cyl}}{\left(\displaystyle \frac{z}{H}\right)}^{-2\beta },\end{eqnarray} \tag{ 19 }$

where the vertical magnetic field is assumed to be constant for a given z and the wave front of the Alfvén wave is assumed to be perpendicular to the vertical axis for simplicity. The time for which the Alfvén wave reaches z is calculated to be

$\begin{eqnarray}&&t(z)={\int }_{0}^{z}\displaystyle \frac{{dz}}{{v}_{{\rm{A}}}(z)},\end{eqnarray} \tag{ 20 }$

where the interval of integration is approximated as [0 : z] from the exact formula of [H_cyl : z]. Note that this simplification is also implicit in Equation (11). The integration yields a solution for z as a function of time as

$\begin{eqnarray}&&z(t)=\displaystyle \frac{(1+2\beta ){H}^{2\beta }{B}_{\mathrm{cyl}}}{{(4\pi {\rho }_{\mathrm{ext}})}^{1/2}}t.\end{eqnarray} \tag{ 21 }$

Then the moment of inertia of the outer region as a function of time is calculated to be

$\begin{eqnarray}&&{I}_{\mathrm{ext}}(t)={\int }_{0}^{t}{\rho }_{\mathrm{ext}}{R}_{\mathrm{ext}}{(z(t))}^{4}{v}_{{\rm{A}}}(z(t)){dt},\end{eqnarray} \tag{ 22 }$

where the interval of integration is again approximated to [0 : t], not to [t₀ : t], where t₀ is defined by z(t₀) = H₀. By using Equation (8), we finally obtain the magnetic braking timescale in a gradually expanding magnetic field configuration as

$\begin{eqnarray}&&{t}_{{\rm{b}},\mathrm{power}}=\displaystyle \frac{1}{1+2\beta }{\left[(1+4\beta )\left(\displaystyle \frac{{\rho }_{\mathrm{cyl}}}{{\rho }_{\mathrm{ext}}}\right)\right]}^{\tfrac{1+2\beta }{1+4\beta }}\displaystyle \frac{{H}_{\mathrm{cyl}}}{{v}_{{\rm{A}}}({H}_{\mathrm{cyl}})}.\end{eqnarray} \tag{ 23 }$

This can be directly compared with Equations (9) and (11). Using the mass-to-flux ratio, Equation (23) is expressed also as

$\begin{eqnarray}&&{t}_{{\rm{b}},\mathrm{power}}=\displaystyle \frac{{(1+4\beta )}^{\tfrac{1+2\beta }{1+4\beta }}}{1+2\beta }{\left(\displaystyle \frac{{\rho }_{\mathrm{cyl}}}{{\rho }_{\mathrm{ext}}}\right)}^{\tfrac{-2\beta }{1+4\beta }}{\left(\displaystyle \frac{\pi }{{\rho }_{\mathrm{ext}}}\right)}^{1/2}\displaystyle \frac{M}{{\rm{\Phi }}}.\end{eqnarray} \tag{ 24 }$

When β = 0, Equations (23) and (24) reduce to Equations (9) and (10), respectively. Also, t_b,power satisfies t_b,power < t_b,∥ because (ρ_cyl/ρ_ext)^{−2β/(1+4β)} < 1. The magnetic braking timescale becomes shorter because of the expansion of the flux tube.

The ratio of t_b,perp to t_b,power is given by

$\begin{eqnarray}&&\displaystyle \frac{{t}_{{\rm{b}},\perp }}{{t}_{{\rm{b}},\mathrm{power}}}=\displaystyle \frac{2+4\beta }{{(1+4\beta )}^{\tfrac{1+2\beta }{1+4\beta }}}{\left(\displaystyle \frac{{\rho }_{\mathrm{ext}}}{{\rho }_{\mathrm{cyl}}}\right)}^{\tfrac{1}{2+8\beta }}.\end{eqnarray} \tag{ 25 }$

This indicates that t_b,power is not always smaller than t_b,⊥ and that whether ${t}_{{\rm{b}},\perp }/{t}_{{\rm{b}},\mathrm{power}}$ is smaller or greater than unity depends on β and on the density ratio. As β increases, the contribution of the density ratio diminishes quickly owing to the strong dependence of the power-law index on β in the term ${\left(\tfrac{{\rho }_{\mathrm{cyl}}}{{\rho }_{\mathrm{ext}}}\right)}^{\tfrac{1}{2+8\beta }}$ . The numerical factor $\tfrac{2+4\beta }{{(1+4\beta )}^{\tfrac{1+2\beta }{1+4\beta }}}$ is greater than unity for β > 0, and t_b,⊥ can become greater than t_b,power. However, the numerical factor depends on our assumptions and has a large uncertainty. Therefore, we cannot directly conclude whether the magnetic braking timescale of parallel configuration is longer or shorter than that of perpendicular configuration on the basis of the analytic estimate.

Although the above analytic discussion has provided physical insight into magnetic braking, the modeling has involved several assumptions and simplifications that can affect even the qualitative results. Multidimensional simulations are therefore essential to achieve a conclusive understanding of the effect of magnetic braking in collapsing cloud cores.

Does Misalignment between Magnetic Field and Angular Momentum Enhance or Suppress Circumstellar Disk Formation?

Article metrics

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Numerical Method and Initial Conditions

2.1. Numerical Method

2.2. Initial and Boundary Conditions

3. Results

3.1. Density and Angular Momentum Structures

3.1.1. Ideal MHD Simulations

3.1.2. Non-ideal MHD Simulations

3.2. Evolution of the Angular Momentum of the Central Region

3.2.1. Ideal MHD Simulations

3.2.2. Non-ideal MHD Simulations

3.2.3. Evolution of the Central Density

3.2.4. Summary of Section 3.2

3.3. Angular Momentum Evolution of the Fluid Element in the Central Region

3.3.1. Initial Locations of the Gas Falling into the Central Region

3.3.2. Evolution of the Angular Momentum of the Fluid Elements in the Central Region

3.3.3. Evolution of the Angular Momentum of the Fluid Elements Initially in the Spherical Region

4. Discussion

4.1. Three Mechanisms for the Evolution of the Angular Momentum during the Cloud Core Collapse and Their Characteristics

4.2. Importance of the Non-ideal Effect

4.3. Comparison with Previous Studies

5. Conclusions

Appendix: Analytical Discussion of the Magnetic Braking

Does Misalignment between Magnetic Field and Angular Momentum Enhance or Suppress Circumstellar Disk Formation?

Article metrics

Share this article

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Numerical Method and Initial Conditions

2.1. Numerical Method

2.2. Initial and Boundary Conditions

3. Results

3.1. Density and Angular Momentum Structures

3.1.1. Ideal MHD Simulations

3.1.2. Non-ideal MHD Simulations

3.2. Evolution of the Angular Momentum of the Central Region

3.2.1. Ideal MHD Simulations

3.2.2. Non-ideal MHD Simulations

3.2.3. Evolution of the Central Density

3.2.4. Summary of Section 3.2

3.3. Angular Momentum Evolution of the Fluid Element in the Central Region

3.3.1. Initial Locations of the Gas Falling into the Central Region

3.3.2. Evolution of the Angular Momentum of the Fluid Elements in the Central Region

3.3.3. Evolution of the Angular Momentum of the Fluid Elements Initially in the Spherical Region

4. Discussion

4.1. Three Mechanisms for the Evolution of the Angular Momentum during the Cloud Core Collapse and Their Characteristics

4.2. Importance of the Non-ideal Effect

4.3. Comparison with Previous Studies

5. Conclusions

Appendix: Analytical Discussion of the Magnetic Braking