Does misalignment between magnetic field and angular momentum enhance or suppress circumstellar disk formation?

The effect of misalignment between the magnetic field $\magB$ and the angular momentum $\Jang$ of molecular cloud cores on the angular momentum evolution during the gravitational collapse is investigated by ideal and non-ideal MHD simulations. For the non-ideal effect, we consider the ohmic and ambipolar diffusion. Previous studies that considered the misalignment reported qualitatively contradicting results. Magnetic braking was reported as being either strengthened or weakened by misalignment in different studies. We conducted simulations of cloud-core collapse by varying the stability parameter $\alpha$ (the ratio of the thermal to gravitational energy of the core) with and without including magnetic diffusion The non-ideal MHD simulations show the central angular momentum of the core with $\theta=0^\circ$ ($\Jang \parallel \magB$) being always greater than that with $\theta=90^\circ$ ($\Jang \perp \magB$), independently of $\alpha$, meaning that circumstellar disks form more easily form in a core with $\theta=0^\circ$. The ideal MHD simulations, in contrast, show the the central angular momentum of the core with $\theta=90^\circ$ being greater than with $\theta=0^\circ$ for small $\alpha$, and is smaller for large $\alpha$. Inspection of the angular momentum evolution of the fluid elements reveals three mechanisms contributing to the evolution of the angular momentum: (i) magnetic braking in the isothermal collapse phase, (ii) selective accretion of the rapidly (for $\theta=90^\circ$ ) or slowly (for $\theta=0^\circ$) rotating fluid elements to the central region, and (iii) magnetic braking in the first-core and the disk. The difference between the ideal and non-ideal simulations arises from the different efficiencies of (iii).


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, 1995Allen et al. 2003;Mellon & Li 2008;Inutsuka et al. 2010;Machida et al. 2011;Dapp et al. 2012;Braiding & Wardle 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 (Price & Bate 2007;Machida et al. 2007;Mellon & Li 2008;Dapp & Basu 2010;Li et al. 2011;Dapp et al. 2012;Tomida et al. 2015;Machida et al. 2014;Lewis et al. 2015;Tsukamoto et al. 2015a; Lewis & Bate 2017;Tsukamoto et al. 2017a). 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 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 three-dimensional simulations of the gravitational collapse of a molecular cloud core with B ∦ Jang (where Jang and 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). 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 (ρ 10 −13 g 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) suggests that the magnetic diffusion eliminates the effect of misalignment on the magnetic-braking efficiency. Their results also imply that the difference between the magneticbraking efficiency in the parallel and perpendicular configurations observed in ideal MHD simulations originates in the high-density region ρ 10 −13 g cm −3 because magnetic diffusion becomes dynamically important only in this region (Umebayashi & Nakano 1990;Nishi et al. 1991;Nakano et al. 2002;Tsukamoto et al. 2015b;Tomida et al. 2015). 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 • was not clear. Wurster et al. (2017) reported that the separation of the binary formed in the core with the perpendicular configura-tion 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;Davidson et al. 2014;Hull et al. 2014;Segura-Cox et al. 2015;Ward-Thompson et al. 2017;Hull et al. 2017b,a;Alves et al. 2018;Soam et al. 2018;Maury et al. 2018;Galametz et al. 2018;Cox et al. 2018;Sadavoy et al. 2018;Kwon et al. 2018a,b).
Furthermore, recent high-resolution observations of Class 0/I YSOs have reported the existence of Keplerian disk with size of several 10s to 100s AU scales (Tobin et al. 2012;Takakuwa et al. 2012;Murillo et al. 2013;Chou et al. 2014;Ohashi et al. 2014;Aso et al. 2015;Yen et al. 2015Yen et al. , 2017Aso et al. 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 α ≡ E therm /Egrav, where E therm and Egrav 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 §2. The simulation results are then shown in §3. We summarize our numerical results and discuss their implications in §4. Finally, we present our conclusions on the impact of misalignment in §5. Appendix A, reviews the analytic discussion of the magnetic-braking timescale. The derived timescales are used for discussions.

Numerical Method
The simulations solve non-ideal radiation magnetohydrodynamics equations with self-gravity: where ρ is the gas density, P is the gas pressure, B is the magnetic field andB ≡ B/|B|, ηO and ηA are, respectively, the resistivity for ohmic and ambipolar diffusion, Er is the radiation energy, Fr is the radiation flux, Pr is the radiation pressure, Tg is the gas temperature, κP is the Plank mean opacity, e = ρu + 1 2 (ρv 2 + B 2 ) is the total energy (where u is specific internal energy), Φ is the gravitational potential, and ar 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);Wurster et al. (2018a).
To close the equations for radiation transfer, we employed the flux-limited diffusion (FLD) approximation (Levermore & Pomraning 1981), where κR is the Rosseland mean opacity. We use the smoothed particle hydrodynamics (SPH) method (Lucy 1977;Gingold & Monaghan 1977;Monaghan & Lattanzio 1985), the numerical code was developed by the present authors and employed previously (e.g., Tsukamoto & Machida 2011Tsukamoto et al. 2013bTsukamoto et al. , 2015aYoneda 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 & Inut-suka 2013). Radiative transfer is solved implicitly using the method of Whitehouse & Bate (2004); 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 (STS) (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 H2, H, H + , He, He + , 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 H + 3 , HCO + , Mg + , He + , C + , 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 −2 µm and ρ d = 2 g 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 non-thermal ionization by the cosmic ray and the thermal ionization of potassium. The cosmic-ray ionization rate was fixed to be ξCR = 10 −17 s −1 .

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 Erot/Egrav = 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 µ = (M/Φ)/(M/Φ)crit = 4 where Φ = πR 2 core B0 (Rcore is the core radius and B0 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 In the following simulations, we scanned over the parameter space comprising two parameters: the ratio α of the thermal energy to the gravitational energy (α ≡ E therm /Egrav) 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 . Thus, the rotation energy is also different between cores with different value 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 4 for α = 0.2, 0.4, 0.6, θ = 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 × 10 6 SPH particles. We conducted the simulations until the epoch just after protostar formation (when the central density ρc becomes ∼ 10 −2 g cm −3 ). A boundary condition was imposed at Rout = 0.95Rcore and the particles with r > Rout 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 ρ < 4.0 × 10 −17 g cm −3 .

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 non-zero θ, 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 θ = 0 • . 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 sit 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 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 vr ≡ (yvy + zvz)/ 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 Jyz = (J 2 y +J 2 z ) 1/2 . Jyz 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).
In the cores with θ = 45 • (panels (d), (e), and (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.

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 ρ 10 −13 g cm −3 (Umebayashi & Nakano 1990;Nakano et al. 2002;Tomida et al. 2015;Tsukamoto et al. 2015b).
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 after 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 figure 4 indicate that the morphology of the angular-momentum map is similar for the ideal and non-ideal MHD simulations although strength differ On the other hand, the outflow morphology is different. For example, the panels (d) and (e) of figures 2 and figure 4 show that the gas in the central region of |y|, |z| < 100AU in nonideal 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). Table 1. Model names and simulation parameters: α = E therm /Egrav the initial radius R 0 , the initial density ρ 0 , the initial angular velocity Ω 0 , the initial magnetic field strength B 0 , the angle between the initial magnetic field and the initial angular momentum vector of cloud cores θ. The last column indicates whether the ohmic and ambipolar diffusions are included ("Yes") or not ("No"). Figure 5 shows the evolution of the mean angular momentum of the central region (the region with ρ > 10 −12 g 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 ρc ∼ 10 −2 g 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).

Ideal MHD simulations
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 Ciardi (2009) andJoos 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 g cm −3 ρc 10 −10 g 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 g cm −3 ρc 10 −8 g cm −3 . The decrease of the central angular momen-tum 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 because 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. , 2016. Figure 6 shows the evolution of the mean angular momentum of the central region with ρ > 10 −12 g cm −3 in the nonideal 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 α.

Non-Ideal MHD simulations
Because the magnetic diffusions are effective in the relatively high density of ρ 10 −13 g 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 ρc > 10 −10 g 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 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 ρc ∼ 10 −2 g cm −3 ) in the simulations of the ideal MHD models (Model I1-I9). The black and red lines represent the contours of the density and of vr ≡ (yvy +zvz)/ y 2 + z 2 = 0, respectively. Black contour levels are ρ = 10 −17 , 10 −16 , 10 −15 , 10 −14 , 10 −13 , and 10 −12 g 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. θ = 45 • are small. For all values of α, the central angular momenta differ between the models with θ = 0 • and θ = 45 • are 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 θ, and consistent with the results of ideal MHD simulations. This suggests that the difference in α mainly affects magnetic braking in the isothermal-collapse phase (ρc < 10 −13 g cm −3 ) in which the magnetic diffusion does not play a role.

Evolution of the central density
In figure 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 similar evolutionary stage, for example, at the epochs of the first-and (ρc ∼ 10 −13 g cm −3 ) and second-core (or protostar) formation (ρc ∼ 10 −2 g 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 y + J 2 z ) 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 ρc ∼ 10 −2 g cm −3 ) in the simulations of the ideal MHD models (Model I1-I9). The black and red lines represent contours of Jyz of vr ≡ (yvy + zvz)/ y 2 + z 2 = 0, respectively. The region enclosed by the red contours is outflow. The black contour levels are Jyz = 10 4 , 10 5 , and 10 6 g cm −1 s −1 . The red arrows show the velocity field. free-fall time by a factor between 1.1 and 1.4. where the free-fall time is t ff = 1.6 × 10 4 yr for α = 0.2, t ff = 4.6 × 10 4 yr for α = 0.4, and t ff = 8.5 × 10 4 yr for α = 0.6. The figure shows that the epochs of the first-(ρc ∼ 10 −13 g cm −3 ) and second-core formations (ρc ∼ 10 −2 g cm −3 ) depend on θ, even for the same α. In the ideal MHD simulations, the models with a greater θ require a longer time for the firstand second-core formations. For example, in the ideal MHD simulations with α = 0.6, the second-core formation epoch differs by ∼ 5 × 10 3 yr between the simulations with θ = 0 • and θ = 90 • . This introduces complexity when we compare the simulations that have different θ value even when their initial cores had the same α and βrot. 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 ρc 10 −14 g cm −3 between the ideal and non-ideal MHD simulations.
The epoch of the second-core formation in the nonideal 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 • with α = 0.2 is due to 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 rcent is estimated to be rcent ∼ j 2 /(GM ) ∼ 4.5(j/(3 × 10 19 cm 2 s −1 )) 2 (M/(0.1 M )) AU. The centrifugal force with a specific angular momentum of the order of j ∼ 10 19 cm 2 s −1 provides the rotational support against self-gravity at a radius of several AU, which corresponds . 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 ρc ∼ 10 −2 g cm −3 ) in the simulations of the non-ideal MHD models (Model NI1-NI9). The black and red lines represent contours of the density and of vr ≡ (yvy + zvz)/ y 2 + z 2 = 0, respectively. The black contour levels are ρ = 10 −17 , 10 −16 , 10 −15 , 10 −14 , 10 −13 , and 10 −12 g 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.
to the radius of the first-core. Figures 6 and 7 show that the duration of the first-core epoch (10 −13 g cm −3 < ρc < 10 −8 g cm −3 ) is considerably longer in the models with j 10 19 cm 2 s −1 than those with j 10 19 cm 2 s −1 . The difference in the second-collapse epoch arises in consequence.
It is worth noting that the order reversal of the secondcollapse 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.

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 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 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 which employed a large value of α (Matsumoto & Tomisaka 2004) and a small α (Hennebelle & Ciardi 2009;Joos et al. 2012). the results becomes simpler once magnetic diffusion has been incorporated into the simulation. The central angular momentum decreases as θ increases, regardless of the value of α (figure 6). . The black and red lines represent contours of the density and vr ≡ (yvy + zvz)/ y 2 + z 2 = 0, respectively. The region enclosed by the red contours is the outflow. The black contour levels are Jyz = 10 4 , 10 5 , and 10 6 g cm −1 s −1 . The red arrows represent the velocity field. Figure 5. The evolution of the mean specific angular momentum in the region with ρ > 10 −12 g cm −3 as a function of the central density in the ideal MHD simulations (Model I1-I9). The horizontal axis represents the central density ρc which is an increasing function of time (see figure 7). The vertical axis represents the mean specific angular momentum of the region with ρ > 10 −12 g cm −3 , where the mean specific angular momentum is calculated asj(ρ) = ρ>10 −12 g cm −3 ρr × vdV/ ρ>10 −12 g cm −3 ρdV.

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 which determines the "disk" (Masson et al. 2016). However, these kinds of analyses make it is 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 new-born 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 ρc ∼ 10 −2 g cm −3 ). Following the criterion used in the previous section, we define the central region as the region where ρ > 10 −12 g cm −3 . At the end of each simulation run, 10 4 SPH particles were sampled from the central region. We then traced back their trajectories backward and determined their initial positions. 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 Lorenz 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 2 + y 2 ), 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 which initially have a large angular momentum selectively accreted onto the central region when θ = 90 • . If we only investigate the central angular momentum, this non-spherical 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., (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 α. 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 figure  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 (where indicated with arrows (i) in figure 8)) is the non-spherical accretion discussed in the previous subsection. Figure 9 indicates that the mean specific angular momentum decreases in the two different density ranges. A decrease occurs in ρ 10 −13 g 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.

Evolution of the angular momentum of the fluid elements in the central region
The other decrease occurs in ρ 10 −11 g 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 two-dimensional 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 free-fall 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 inves-tigating 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 4 SPH particles were sampled from the region with ρ > 10 −12 g 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 (ρ 10 −13 g 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 ρ 10 −13 g 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 ρ 10 −13 g 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 is also weak. As a result, the efficiency of magnetic braking is reduced in the high-density region compared to the ideal MHD simulations.

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 which is selected from the initial cloud core. Let the sphere be parameterized by the enclosed mass, We chose a sphere with Me(r) = 5 × 10 −2 M , i.e., the innermost region of the entire cloud core. Then 10 4 SPH Figure 8. Initial distribution of the SPH particles that have fallen into the central region (with density greater than 10 −12 g cm −3 at the end of the simulation), obtained from the ideal MHD simulations (Model 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 2 + y 2 ).  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.
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 θ value is small. Therefore, in the simulations with a 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 4 ) is plotted with dashed-dotted lines. The lines show that the angular momentum is almost constant during the isothermalcollapse phase. This confirms that the decrease of the angular momentum observed in MHD simulations is indeed caused by the magnetic field.

Three mechanisms for the evolution of the angular momentum during the cloud core collapse and their characteristics
In §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; (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 which 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 A1 is smaller than that with θ = 90 • and a field geometry of panel (b) of figure A1 by a factor of (Rc/Rcore)( 1) times, where Rc and Rcore 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 Appendix A, 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 isothermalcollapse phase decreases as α decreases (figures 9 and 11). In this study, we parameterized the initial cloud core with α ≡ E therm /Egrav. Because the initial temperature is fixed to 10 K, a small α corresponds to a large gravitational energy and short free-fall 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 becomes 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 5 years 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 highdensity 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 free-fall time of the initial cloud core is not the characteristic time-scale. Owing 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 hour-glass like magnetic field geometry as indicated by equations (A16) and (A18). This, along with non-spherical 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.

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 diffusion are effective in the region with ρ 10 −13 g cm −3 and in ρ 10 −11 g cm −3 , respectively (Dapp et al. 2012;Tsukamoto et al. 2015b;Tomida et al. 2015;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 nonideal 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 nonideal 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 to borne in mind when comparing the previous studies with non-ideal effects (e.g., Tomida et al. 2015;Tsukamoto et al. 2015b;Wurster et al. 2016;Masson et al. 2016;Zhao et al. 2016;Wurster et al. 2018a;Vaytet et al. 2018;Zhao et al. 2018b) 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; Wurster et al. 2018a,b). 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 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., ρ > 10 −12 g cm −3 ) of the core with θ = 180 • is an order of magnitude greater 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 value 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.

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. Ciardi (2009) andJoos et al. (2012) reported that magnetic braking is more efficient in a core with θ = 0 • than with θ = 90 • , in apparent contradiction with by Matsumoto & Tomisaka (2004). Note, however, that the difference in results of Hennebelle & Ciardi (2009), Joos et al. (2012 and Matsumoto & Tomisaka (2004) not contradictory but originate from different samplings of the parameter space, as we discuss bellow. Hennebelle & Ciardi (2009), 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 mentions of the configura-tion of the magnetic field was given in Joos et al. (2012). Their cores have α = 0.25. Using a barotropic equation of state 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 nc > 10 10 cm −3 , and that the mean angular momentum, in the range including the low density region (nc > 10 8 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 significant fraction of the angular momentum is removed during the isothermalcollapse 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, Ciardi (2009) andJoos 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 hour-glass like magnetic field (see Appendix A), the central angular momentum in the core with θ = 0 • is smaller than that with θ = 90 • . The selective accretion may also assist that 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 B and Jang (θ < 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 equation of state 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 equation of state, or numerical schemes. We think the analysis of the evolution of the angular momentum evolution 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 element as described in this paper.

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.

ACKNOWLEDGMENTS
We thank Tomoyuki Hanawa and Tomoaki Matsumoto for fruitful discussion. We also thank anonymous referee for helpful comments which greatly improve 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) to us. The computations were performed on a parallel computer, XC40/XC50 system at Center for Computational Astrophysics of National Astronomical observatory of Japan. This work is supported by JSPS KAKENHI grant number 17KK0096, 17K05387, 17H06360, 17H02869, 18K13581.

APPENDIX A: ANALYTICAL DISCUSSION OF THE MAGNETIC BRAKING
In this section, we derive the magnetic-braking timescales based on the analytical model of Mouschovias & Paleologou (1979 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 configuration. Mouschovias & Paleologou (1979 showed that the magnetic-braking timescale t b of a central region with a moment of inertia Ic 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 Iext(t b ) equals to Ic. This condition is given by This then determine 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 (Jang B or θ = 0 • ) can be calculated as follows. The central collapsing region is modeled as a uniform cylinder with a density ρ cyl , radius R cyl , and height 2H 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) in figure A1). The moments of inertia of the central cylinder and outer envelope are given by Ic = πρ cyl R 4 cyl H cyl and Iext(t b ) = πρextR 4 cyl vAt b , respectively, where vA denotes the Alfvén velocity in the outer envelope and vAt b corresponds to the height of the swept outer envelope. By solving equation (A1), t b is calculated to be (Mouschovias 1985). Using the mass M = 2πρ cyl R 2 cyl H cyl and the magnetic flux Φ = πR 2 cyl B of the cylinder, we can reduce t b, to

Does misalignment between magnetic field and angular momentum enhance or suppress circumstellar disk formatio
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.
Magnetic-braking timescale in a simple geometry with a perpendicular configuration (Jang ⊥ 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 ρ cyl , R cyl , and 2H cyl . The magnetic field is assumed to be perpendicular to the rotation axis and axisymmetrical (panel (b) in figure A1) where B(r) ∝ r −1 . The density of the outer envelope ρext is again assumed to be constant (panel (b) of figure A1).
The moments of inertia of the central cylinder and outer envelope are given by Ic = 4πρ cyl H cyl R 4 cyl and Iext = 4πρextH cyl (R 4 ext − R 4 cyl ). Using the condition Ic = Iext, the magnetic-braking timescale for the perpendicular configuration is calculated to be (Mouschovias 1985), where we assumed ρ cyl /ρext 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 vA ∝ 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 Thus the timescale for perpendicular configuration is shorter than that for the parallel one because ρ cyl ρext and the magnetic braking in the former is stronger than in the later. 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 hour-glasslike shape as a result of radial dragging according to the theoretical work of (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 Rext just above the central cylinder (panel (c) in figure A1). If we ignore the moment of inertia of the transition layer, Iext(t b ) is given by Using Ic = πρ cyl R 4 cyl H cyl and equation (A7), we obtain the magnetic-braking timescale of the disk with an hour-glass magnetic field geometry (Mouschovias 1985) as where is assumed because of the conservation of magnetic flux. According to equation (A8), the magnetic-braking timescale can be shorter than t b, in equation (A3) because of the factor (R cyl /Rext) 2 (< 1). Joos et al. (2012) argued that the magnetic-braking timescale becomes smaller in the parallel than perpendicular configuration by virtue of equation (A8). They assumed Rext = Rcore, where Rcore is the initial core radius, and ρ(r) ∝ r −2 . Then, assuming that the mean external density and the central density can be approximated in terms of the volume-averaged means as, respectively, ρ cyl ∼ ρ(R cyl ) and ρext ∼ ρ(Rext), the ratio of the timescales is estimated to be Using this relation, they concluded that the braking timescale in flared parallel configuration is shorter than in the perpendicular configuration because R cyl < Rcore.
However, it is uncertain whether the assumption of an infinitely thin magnetic fluxtube 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 /Rext) 2 . This assumption corresponds to discontinuous decrease of poloidal magnetic field strength (equation (A9)) 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 Rext.
We should also note that validity of the assumption of Rext = Rcore is unclear, because Rext is the radius of the external flux tube and is not related to that of the core Rcore, 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). Rext can be Rext Rcore even at the vertical height of z = Rcore 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 (A8) 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 which expand in a power-law fashion. Our result will show that whether the timescale becomes shorter or not