MAGNETIC BRAKING AND PROTOSTELLAR DISK FORMATION: AMBIPOLAR DIFFUSION

Disk formation is an integral part of star formation that has been studied for a long time. Early works concentrated on the collapse of rotating, nonmagnetic cores (e.g., Terebey et al. 1984). Bodenheimer (1995) reviewed these works, and noted a number of unresolved problems (see also Boss 1998). Topping the list is the effect of magnetic braking on disk formation.

The magnetic properties of star-forming dense cores are reasonably well constrained. There is now ample evidence for ordered magnetic fields on the core scale from polarization of dust emission (e.g., Ward-Thompson et al. 2000). A spectacular recent example is the Submillimeter Array (SMA) observation of NGC 1333 IRS 4A, which shows a pinched magnetic configuration on the 10³ AU scale (Girart et al. 2006). The magnitude of the magnetic field is harder to determine. Troland & Crutcher (2008) carried out an extensive Zeeman survey of dark cloud cores in OH, and determined a mean mass-to-flux ratio of λ ∼ 4.8 in units of the critical value (2πG^1/2)⁻¹ (Nakano & Nakamura 1978); correction for uncertain projection effects may bring this value to ∼2. The implication is that dense cores are typically moderately strongly magnetized, with a dimensionless mass-to-flux ratio of a few to several. A relatively low value of λ ∼ 2–3 is expected for dense cores formed out of magnetically subcritical clouds (with λ < 1) through ambipolar diffusion (e.g., Lizano & Shu 1989; Basu & Mouschovias 1994; Nakamura & Li 2008; Adams & Shu 2007). The values of λ are expected to have a wider spread for cores formed out of turbulent compression. Nevertheless, the cores tend to be more magnetized relative to their masses than the cloud as a whole, because only a fraction of the mass along a flux tube that threads the cloud is compressed into the core (e.g., Dib et al. 2007; Tilley & Pudritz 2005). It is unlikely for the cores to have a value of λ more than several, unless the cloud as a whole is magnetized to an unrealistically low level.

There have been a number of recent calculations that included magnetic fields in the collapse of rotating cores, focusing mostly on the early phase before the mass of the central object reaches the stellar range (e.g., Tomisaka 1998; Machida et al. 2006; Banerjee & Pudritz 2006; Price & Bate 2007). Other studies have concentrated on the angular momentum evolution and disk formation in the main accretion phase (Allen et al. 2003; Galli et al. 2006; Mellon & Li 2008, hereafter Paper I; Hennebelle & Fromang 2008). These studies find that magnetic braking efficiently removes angular momentum from the collapsing matter, preventing a large, 10² AU scale rotationally supported disk from forming in realistic cores in the ideal MHD limit. However, such disks are routinely observed around young stars, at least in relatively late, Class I and II phases. How can the disk be saved?

Molecular cloud cores are lightly ionized, so perfect coupling between the magnetic field and matter assumed in ideal MHD is not expected. As the density increases, the coupling is weakened first by ambipolar diffusion, then the Hall effect, then Ohmic dissipation. In this paper, we will concentrate on ambipolar diffusion as a first step. The effect of ambipolar diffusion on magnetic braking in the formation (Basu & Mouschovias 1994) and collapse (Krasnopolsky & Königl 2002) of rotating cores has been examined previously using the thin disk approximation. Krasnopolsky & Königl (2002) parameterized the braking strength and found that disk formation can be suppressed in the ambipolar diffusion limit when the braking parameter is large. We use two-dimensional models which allow direct calculation of the strength of the magnetic braking. We find that ambipolar diffusion does not weaken the magnetic braking enough to allow rotationally supported disks to form under realistic conditions. The implication is that additional processes must be found to save the disk; two possibilities are discussed in Section 5.

We study disk formation and evolution during the main accretion phase after a central object has formed. The problem setup is identical to that of Paper I (to which we refer the reader for details), except for the inclusion of ambipolar diffusion. As in Allen et al. (2003), we idealize the initial configuration of the main accretion phase as a rotating self-similar singular isothermal toroid supported against self-gravity partially by thermal pressure and partially by magnetic fields (see Figure 1 of Paper I for an example). This initial configuration represents the pivotal state separating the prestellar and protostellar phases, which is reached via ambipolar diffusion initiated formation of dense cores. The flattening of this pivotal state is due to anisotropic magnetic support; it depends on the degree of magnetization. This pivotal state is in force balance with no initial radial velocity. Tests performed in Paper I found that an initial subsonic infall does not change the results qualitatively. The toroid is characterized by three parameters, the isothermal sound speed (c_s), mass-to-flux ratio (λ), and rotational speed (v₀), which is taken to be spatially constant to preserve the self-similarity of the configuration.

Our treatment of ambipolar diffusion follows that of Shu (1991), where the magnetic field is assumed to be tied to ions, which drift relative to neutrals with a velocity

$$\begin{equation} {\bf v}_d={\bf v}_i-{\bf v}_n=\frac{1}{4\pi \gamma \rho _n\rho _i} (\nabla \times {\bf B}) \times {\bf B}, \end{equation} \tag{ 1 }$$

where γ = 3.5 × 10¹³ cm³ g⁻¹ s⁻¹ is the ion–neutral drag coefficient and ρ_n (ρ_i) and v_n (v_i) are the density and velocity of the neutrals (ions), respectively. We adopt an ion density of the form

$$\begin{equation} \rho _i=C \rho _n^{1/2} \end{equation} \tag{ 2 }$$

with the coefficient

$$\begin{equation} C=3\times 10^{-16}\left(\frac{\zeta }{10^{-17} {\rm s}^{-1}} \right)^{1/2}\mbox{cm}^{-3/2}\, {\rm g}^{1/2}, \end{equation} \tag{ 3 }$$

where ζ is the cosmic ray ionization rate (Elmegreen 1979). This simple form has the advantage of preserving the self-similarity of the collapse, which provides a powerful check on our numerically obtained solutions. We implement the ambipolar diffusion into the ZeusMP MHD code (Hayes et al. 2006) using the fully explicit method of Mac Low et al. (1995). As is well known, for such explicit methods to be stable, the limiting time step must be proportional to the grid size squared, which puts stringent constraints on spatial resolution.

We solve the governing MHD equations including ambipolar diffusion in a spherical coordinate system (r, θ, ϕ) assuming axisymmetry. The grid spacing is logarithmic in radius r and constant in polar angle θ, with a computational domain extending radially from 10¹⁴ to 2 × 10¹⁷ cm and 0 to π in angle. There are 120 grid points in the radial direction, and 60 angular points. The smallest grid size is ∼10¹³ cm, near the inner boundary. The standard outflow condition is imposed at the outer radial boundary. The inner boundary is a modified outflow boundary with the mass accreted across the boundary added to the central point mass and a torque free condition imposed on the magnetic field (i.e., B_ϕ = 0, see Paper I). We adopt the oft-used broken power-law equation of state that is isothermal below ρ = 10⁻¹³ g cm⁻³ and adiabatic with γ = 7/5 above.

We first illustrate the effects of ambipolar diffusion on magnetic braking and disk formation using a model with a fiducial cosmic ray ionization rate ζ = 10⁻¹⁷ s⁻¹ and a mass-to-flux ratio λ = 4; the latter is close to the mean value inferred by Troland & Crutcher (2008) in their Arecibo OH Zeeman observations of dark cloud cores.

Other combinations of ζ and λ are considered in the next section. For all models, we set the sound speed to c_s = 0.3 km s⁻¹ and the rotational speed to v₀ = c_s/2 = 0.15 km s⁻¹, which corresponds to an angular speed of 3 km s⁻¹ pc⁻¹, typical of dense cores (Goodman et al. 1993). Our choice of model parameters is motivated by the properties inferred for L1544, arguably the best studied starless core (e.g., Doty et al. 2005).

A quantity of primary importance to star formation is the mass accretion rate, ${\dot{M}}$. In Figure 1, we show ${\dot{M}}$ as a function of time. In our standard model, ${\dot{M}}$ converges to a constant value, after a short period of adjustment. The initial deviation from the constant is due to the zero infall velocity at t = 0 and the small point mass that we put at the center to induce the collapse. The accretion rate quickly approaches a constant value indicating that the solution has reached the expected self-similar state at late times. The converged value of ${\dot{M}}$ is close to that of Shu (1977) for a singular isothermal sphere. One might naively expect the rotation to retard the collapse significantly, leading to a reduced mass accretion rate onto the central object. Apparently, the magnetic braking is strong enough that rotation is not a significant barrier to mass accretion, even in the presence of ambipolar diffusion. The accretion rate is comparable to that in the ideal MHD case (Paper I), suggesting that the strength of braking is not significantly reduced by ambipolar diffusion. The accretion is highly episodic in the ideal MHD case, but is steady in the presence of ambipolar diffusion.

Zoom In Zoom Out Reset image size

Figure 2 shows a snapshot of the standard model at a representative time of 5.85 × 10¹¹ s. As in the singular isothermal sphere, the collapse occurs inside out, with the bulk of the envelope material beyond the radius of c_st = 2 × 10¹⁶ cm remaining nearly static. Inside this radius, the matter accelerates toward the center. The collapsing flow is deflected toward the equator by pinched field lines, forming a pseudodisk (Galli & Shu 1993). As the flow collapses, it is expected to spin up due to conservation of angular momentum. However, the actual rotation speed tends to decrease as the collapsing material moves inward, especially in the equatorial region, indicating significant loss of angular momentum due to magnetic braking. In the evacuated polar cavities, the angular momentum is removed from the system via a low velocity (∼0.5 km s⁻¹) outflow. This outflow is not unique to the ambipolar diffusion case (see Paper I), and is not the bipolar molecular outflow often observed in the early phase of star formation (e.g., Gueth & Guilloteau 1999), which is most likely driven by a fast jet or wind originated interior to our computational space.

Zoom In Zoom Out Reset image size

To illustrate the braking more quantitatively, we plot the equatorial velocity profiles in Figure 3. As one moves from large radii to the center we observe four distinct regions in the radial velocity profile (top panel). Outside of r = 2 × 10¹⁶ cm the matter remains nearly static. In the region 6 × 10¹⁵ cm ≲r ≲ 2 × 10¹⁶ cm the matter accelerates toward the center as gravity dominates the magnetic forces, since the mass-to-flux ratio of the collapsing core is significantly greater than unity. The collapsing material is decelerated in the region 3 × 10¹⁵ cm ≲r ≲ 6 × 10¹⁵ cm, where the outward magnetic forces exceed the inward gravitational pull. This region does not exist in the ideal MHD limit; its presence is due entirely to ambipolar diffusion, which enables the magnetic field lines that would have been dragged into the central object in the ideal MHD limit to pile up outside the object (Li & McKee 1996; Ciolek & Königl 1998). In this decelerating matter, the ion–neutral drift velocity becomes large because of a large Lorentz force. At r ∼ 3 × 10¹⁵ cm the neutral matter begins accelerating toward the center again. However, the ions remain nearly static at r ∼ 1.5 × 10¹⁵ cm, indicating that the magnetic field lines barely move in that region. The deceleration region has the same physical origin as the hydromagnetic accretion shock induced by ambipolar diffusion (Li & McKee 1996; Ciolek & Königl 1998). However, since the gas outside of the decelerating region is collapsing subsonically, a hydromagnetic accretion front rather than a shock is formed. At r ≲ 10¹⁵ cm the ions (and magnetic field) begin reaccelerating toward the center. The rapid collapse of both ions and neutrals at the smallest distances from the star clearly suggests the absence of a rotationally supported disk.

Zoom In Zoom Out Reset image size

From the lower panel of Figure 3, it is clear that the rotation speed v_ϕ does not increase with decreasing radius, which is opposite of what is expected based on angular momentum conservation. For example, in the region of accelerating collapse (6 × 10¹⁵ cm ≲r ≲ 2 × 10¹⁶ cm), v_ϕ remains constant while the radius decreases by a factor of ∼3. This indicates significant loss of angular momentum. This loss can be understood by examining the ratio of toroidal to poloidal magnetic fields. We find that the magnetic field lines get more twisted as the collapse accelerates inward. When deceleration occurs, compression increases the poloidal and toroidal field strength, which leads to a stronger braking that is responsible for the steep drop in rotational speed in the deceleration region (the shaded region in Figure 3). By the time the fluid reaccelerates toward the center, there is little angular momentum left in the gas to maintain a significant twist in the magnetic field, resulting in near radial infall to the center. The supersonic infall velocity and small rotation speed at small radii leaves no doubt that rotationally supported disk formation is suppressed in our standard model.

Our collapse solution depends on the cosmic ray ionization rate, which is somewhat uncertain. For dense cores, the typical values lie in a range between (1–5)× 10⁻¹⁷ s⁻¹ (e.g., Maret & Bergin 2007), although a wider spread is also possible (Caselli et al. 1998). To examine the sensitivity of our results to the ionization rate, we vary ζ by a factor of 10 in either direction.

The higher ionization rate case has a larger ionization fraction by a factor of $\sqrt{10}$, which leads to stronger coupling between the matter and magnetic field. Its collapse is qualitatively similar to the standard model, with a few quantitative differences. The stronger coupling allows the magnetic fields to be brought closer to the central object, leaving a larger magnetic flux in a smaller region. As a result, the deceleration region is closer to the central object (see the dashed lines of the top panel of Figure 4). The stronger field and stronger coupling keeps the ion radial speed low for a more extended region during the reacceleration. The stronger field also results in enhanced braking which removes essentially all the angular momentum at small radii (see the dashed lines of the second panel, noting that the rotation speeds for ions and neutrals are practically indistinguishable). Increasing the ionization rate further suppresses rotationally supported disk formation.

Zoom In Zoom Out Reset image size

The lower ionization rate yields an ionization fraction smaller by a factor of $\sqrt{10}$, which leads to weaker coupling between the matter and magnetic field. The ions start to decelerate earlier, due to the pile up of magnetic flux at larger radii. However, the neutrals slip through the ions without being appreciably decelerated due to weak coupling, resulting in faster diffusion. This results in a faster collapse of the neutrals and a slower collapse of the ions compared to the standard model (see the dot-dashed lines in the top panel). The increased diffusion requires a smaller ambipolar diffusion time step, which increases computational time.¹ The weak coupling and reduced magnetic field strength decrease the strength of magnetic braking, allowing the neutrals to spin up as they collapse. However, there is still enough braking to prevent a rotationally supported disk from forming inside our computational domain, as evidenced by the rapid radial collapse at small radii. We conclude that for moderately strongly magnetized core of λ = 4, disk formation is suppressed for realistic values of the cosmic ray ionization rate.

The mass-to-flux ratio of dense cores is uncertain. Our standard model is based on the mean mass-to-flux ratio inferred by Troland & Crutcher (2008), but many cores have only lower limits on the mass-to-flux ratio. More highly magnetized cores have a higher degree of flattening, which would reduce the mass of material in the polar regions. The reduction tends to lower the efficiency of magnetic braking. However, the increase in the poloidal field strength will increase the efficiency of magnetic braking. Allen et al. (2003) found stronger braking in more highly magnetized (thus more flattened) cores (with mass-to-flux ratios as low as 2), indicating that the increase in braking efficiency due to increased poloidal field strength is more than enough to offset the reduction due to flattening. We investigate the collapse model in more weakly magnetized cores with λ = 13.3 (same as the standard ideal MHD case of Paper I), which corresponds to a field strength of 7.35 μG on the scale of the typical core radius of 0.05 pc. We take this case as a lower limit of realistic magnetic field strengths based on the median field strength of cold neutral H i structures of ∼6 μG inferred by Heiles & Troland (2005).

We again consider three values for the ionization rate ζ = 10⁻¹⁶, 10⁻¹⁷, and 10⁻¹⁸ s⁻¹, concentrating first on the case with the fiducial value of 10⁻¹⁷ s⁻¹. In contrast to the standard model with the same ζ, there are two equatorial deceleration regions in this weaker field case instead of one (see the solid lines in the third panel of Figure 4). The accelerating infall is first decelerated near r ∼ 10¹⁶ cm. This deceleration is not a new feature caused by ambipolar diffusion; it was already present in our previous ideal MHD calculation: it is a "magnetic barrier" caused by the bunching of magnetic field lines at the interface between the collapsing envelope and an expanding, magnetic braking-driven bubble (see Figure 3 of Paper I). Inside the barrier, the material recollapses inward, at an increasingly high speed, until a radius of ∼10¹⁵ cm, where it is decelerated for a second time. The second deceleration is caused by ambipolar diffusion, which allows magnetic flux to pile up at small radii, instead of being dragged into the center. As in the standard model, the magnetic pressure in the deceleration region becomes comparable to the ram pressure of the collapsing material. However, in this model the deceleration region occurs at a smaller radius, allowing the infalling gas to collapse supersonically, causing a hydromagnetic shock instead of a front. The slowdown of infall compresses the collapsing material, leading to a stronger poloidal and toroidal magnetic field that enhance the rate of magnetic braking. The enhancement is responsible for the dip on the rotational speed profile (see the solid lines in the bottom panel). Interior to the deceleration region, the collapsing flow accelerates for a third time, spinning up as it collapses. Nevertheless, there is insufficient angular momentum left in the material for a rotationally supported disk to form before the inner boundary is reached.

The dynamics of the more ionized, ζ = 10⁻¹⁶ s⁻¹ case is similar to that of the ζ = 10⁻¹⁷ s⁻¹ case down to the magnetic diffusion-induced second deceleration near ∼10¹⁵ cm. Near this region and interior to it, the infall speed is also similar, but the rotational speed is drastically different (see the dashed lines in the bottom two panels). The field compression in the deceleration region increases the rate of magnetic braking to such a degree as to cause a counterrotation in the better coupled case. The better field-matter coupling and stronger magnetic field lead to efficient removal of angular momentum at small radii, suppressing disk formation completely.

Only in the most weakly ionized case of ζ = 10⁻¹⁸ s⁻¹ does a rotationally support disk form. The disk is most obvious in the radial velocity curve (dot-dashed lines in the third panel). As in the more strongly coupled cases, there is accelerating collapse at large radii. The collapse is only weakly slowed by the magnetic barrier, because the neutrals are weakly coupled to the field. The ions and neutrals are reaccelerated, then decelerated for a second time. Unlike the more strongly coupled cases, this second deceleration is not due to the ambipolar diffusion-enabled pile up of magnetic flux, but is caused by centrifugal force from rapid rotation (see the dot-dashed lines in the bottom panel). The combination of a weak magnetic field and weak matter-field coupling renders the braking too inefficient to remove enough angular momentum to suppress disk formation in this case of extreme parameters.

We have studied the collapse of rotating molecular cloud cores magnetized to moderate degrees, concentrating on a dimensionless mass-to-flux ratio λ = 4, as suggested by recent Zeeman observations (Troland & Crutcher 2008). In the ideal MHD limit, we have previously shown that a slight twist of the field lines in such cores is sufficient to remove essentially all of the angular momentum of the collapsing matter and suppress the formation of rotationally supported disks during the protostellar accretion phase completely (Paper I; see also Galli et al. 2006; Hennebelle & Fromang 2008). One may expect ambipolar diffusion to reduce the efficiency of magnetic braking and potentially enable the disks to form, since the twisting of field lines by the rotation of neutrals will be reduced by diffusion. However, for realistic cosmic ray ionization rates, the small amount of field twisting needed for angular momentum removal drives a drift speed between the ions (and the field lines tied to them) and neutrals in the azimuthal direction that is much smaller than the neutral rotation speed (see the lower panel of Figure 3). By itself, the azimuthal slippage of field lines relative to neutrals does not reduce the strength of the toroidal magnetic field (and the braking rate) significantly. The ambipolar diffusion is expected to modify the poloidal magnetic field as well, especially in regions where the magnetic force is large. The modification is most evident at small radii, where magnetic diffusion has allowed the magnetic flux that would have been dragged onto the central object to occupy an extended region (Li & McKee 1996; Ciolek & Königl 1998; Tassis & Mouschovias 2007). The increased poloidal field strength tends to make the magnetic braking more efficient (Krasnopolsky & Königl 2002). Given these competing effects in opposite directions, it is not at all obvious whether ambipolar diffusion would increase or decrease the overall efficiency of magnetic braking. Our detailed calculations have shown that, for realistic rates of cosmic ray ionization, the overall efficiency is increased, because the azimuthal field slippage is small and the poloidal field trapped at small radii is strong. The efficiency is reduced only for unrealistically low ionization rates, which speed up the azimuthal field slippage and lower the poloidal field strength by enabling the trapped flux to escape to large distances.

There is, however, ample evidence for large, rotationally supported disks on 10² AU scale or more around young stars, at least in Classes I and II, and possibly in Class 0, sources. This begs the question: what is the origin of these disks?

The most obvious possibility is that additional nonideal MHD effects may weaken the magnetic braking further. Ohmic dissipation is expected to be important at high densities (Shu et al. 2006). Nakano et al. (2002) estimates that it dominates other processes at densities above ∼10¹² cm⁻³ (see also Desch & Mouschovias 2001). Since the density in our standard model never exceeds this value, we believe that Ohmic dissipation would not significantly weaken the magnetic braking inside our computation domain (from 10¹⁴ cm <r < 2 × 10¹⁷ cm). However, Ohmic dissipation may still enable a rotationally supported disk to form at a smaller radius where the density is higher, as demonstrated explicitly by Machida et al. (2008) using a resistive MHD code. They found a disk of ∼10¹² cm in radius shortly (∼11 days) after the formation of the so-called second core; the disk was not present in the ideal MHD counterpart. Such Ohmic dissipation-enabled disks can form inside our inner boundary (of 10¹⁴ cm). Since the material crossing the inner boundary of our simulation is typically braked to a rotational speed well below the local Keplerian speed, any disk that may form later should be relatively small in the absence of angular momentum redistribution.

Besides Ohmic dissipation, another nonideal mechanism of magnetic flux diffusion is the Hall effect. In the absence of dust grains, the Hall effect becomes important when the ions begin to decouple from the magnetic field. This occurs when the ion Hall parameter β_i drops below unity (e.g., Nakano 1984). We have checked that β_i > 1 everywhere inside our computational domain for the standard model. We do not expect the Hall effect to change our results significantly, in the absence of dust grains. The inclusion of dust grains would increase the relative importance of the Hall effect compared to ambipolar diffusion at lower densities, though the exact effect depends on the grain size distribution (Wardle & Ng 1999), which is uncertain. In particular, the Hall effect is sensitive to small difference in the abundances of positively and negatively charged grains. It effectively disappears when grains of equal mass and opposite charge dominate the conductivity, which occurs above a density threshold of ∼10¹¹ cm⁻³ according to Wardle & Ng (1999). If the additional nonideal MHD effects do not allow the observed large-scale disks to form, other means of weakening the magnetic braking must be sought.

One possibility is protostellar outflow. Protostellar outflows are observed at all stages of star formation. They are thought to clear away the envelope material as the source ages (Arce & Sargent 2006). As the envelope is the anchor for the magnetic brake in our models, removing it may reduce the strength of braking sufficiently to allow large scale rotationally supported disks to form. This is particularly true for the more evolved Class II and perhaps Class I sources, where little envelope is left to brake the disk. It may not work for Class 0 sources where the outflows tend to be confined to the polar regions (e.g., HH 211, Gueth & Guilloteau 1999). If this is the case, we would not expect these deeply embedded objects to harbor large-scale, rotationally supported disks. They could still have relatively small disks enabled by, for example, Ohmic dissipation. It is difficult to determine the exact amount of rotational support in the equatorial region from current observations of Class 0 sources. The situation will improve when ALMA comes online. An alternative possibility is the absence of a significant envelope in the first place. If, for example, the star-forming core is weakly magnetized but highly flattened along the field lines, the magnetic braking may be dominated by a low density medium external to the core. In such a case, the braking may become too inefficient to suppress disk formation, with or without ambipolar diffusion, as demonstrated explicitly by Krasnopolsky & Königl (2002). How prevalent such cores are in molecular clouds remains to be determined.

We thank Mordicai Mac Low and Jeff Oishi for help in the implementation of ambipolar diffusion in ZEUS MP, and the referee for helpful comments. This work is supported in part by NASA (NNG05GJ49G) and NSF (AST-0307368) grants.