The early history of protostellar disks, outflows, and binary stars
aa r X i v : . [ a s t r o - ph . S R ] N ov The early history of protostellar disks, outflows and binary stars
Dennis F. Duffin
Department of Physics and Astronomy, McMaster University, Hamilton ON L8S 4M1, Canada [email protected] andRalph E. Pudritz
Department of Physics and Astronomy, McMaster University, Hamilton ON L8S 4M1, CanadaOrigins Institute, McMaster University, Hamilton ON L8S 4M1, Canada [email protected]
ABSTRACT
In star formation, magnetic fields act as a cosmic angular momentum extractor thatincreases mass accretion rates onto protostars and in the process, creates spectacularoutflows. However, recently it has been argued that this magnetic brake is so strongthat early protostellar disks – the cradles of planet formation – cannot form. Ourthree-dimensional numerical simulations of the early stages of collapse ( . yr) ofoverdense star–forming clouds form early outflows and have magnetically regulated androtationally dominated disks (inside 10 AU) with high accretion rates, despite the slipof the field through the mostly neutral gas. We find that in three dimensions, magneticfields suppress gravitationally driven instabilities which would otherwise prevent young,well ordered disks from forming. Our simulations have surprising consequences for theearly formation of disks, their density and temperature structure, the mechanism andstructure of early outflows, the flash heating of dust grains through ambipolar diffusion,and the origin of planets and binary stars. Subject headings: binaries: general — circumstellar matter — methods: numerical —stars: formation — stars: general — stars: low-mass, brown dwarfs
1. Introduction
Over the past decade numerical simulations have enabled the exploration of the central physi-cal questions regarding the nature of star formation – from the collapse of an initial dense gaseousmolecular cloud core to the formation of a star and its associated protostellar disk and outflow,and the emergence of a star’s planetary system. The formation of disks and jets during starormation is central to many of these issues, but very little is known about the earliest phasesof their evolution. How is the initial excessive angular momentum associated with the star’snatal core removed – through magnetic braking (Basu & Mouschovias 1994) and then outflows(Banerjee & Pudritz 2006; Pudritz et al. 2007), or by spiral density waves in disks (e.g. Larson2009)? Do multiple stars form through gravitational fragmentation of cores or massive disks(Hosking & Whitworth 2004; Price & Bate 2007; Hennebelle & Teyssier 2008)? What is the sig-nificance of outflows and jets that are launched before most of the mass has collapsed into thedisk (Lynden-Bell 2003; Banerjee & Pudritz 2006) in a wide variety of young stellar systems –from brown dwarfs (Whelan et al. 2005; Machida et al. 2009) to massive stars (Arce et al. 2007;Banerjee & Pudritz 2007)?Extraction of angular momentum by magnetic fields that thread the collapsing gas may betoo efficient, according to recent two-dimensional (Mellon & Li 2008) and three-dimensional ax-isymmetric (Hennebelle & Fromang 2008) simulations, preventing the formation of rotationallydominated disks, even when the effects of imperfect coupling of the field with the gas are included(Mellon & Li 2009). These results seemingly contradict the observations which clearly show thatdisks are present around most if not all young stars, even in environments in which the magneticfield is expected to be strong (Hartmann & Calvet 1995). Hennebelle & Fromang (2008) performedthree-dimensional ideal MHD simulations of collapsing cores using a barotropic equation of state,concluding that no rotationally dominant structure is formed from 10 to 100 AU for highly magne-tized cores. Mellon & Li (2008) performed two-dimensional ideal MHD simulations on collapsingsingular isothermal toroids using a barotropic equation of state and an inner boundary at 6.7 AU(effectively a sink particle), finding that even moderately magnetized disks could not form. Suchtwo-dimensional models impose a high degree of mathematical symmetry and therefore miss theformation of bars and spiral waves in disks.We improve on previous results through this three-dimensional adaptive mesh refinement(AMR) investigation which includes a full treatment of the cooling (Banerjee et al. 2006) and thefinite coupling of the magnetic field to the pre–stellar gas (ambipolar diffusion, Duffin & Pudritz2008). We find that ordered disk–like structures can emerge on scales .
10 AU at early times ( . yr) in magnetized systems. We have omitted a sink particle in this study as they are expected toaffect the solution to within a few sink radii. Without sinks simulations are indeed limited to earlydisks, although they offer a full solution to the region within 10 AU where the heated core forms.
2. Methods
We model an isolated star forming core with a rotating Bonnor–Ebert (BE) sphere (Bonnor1956; Ebert 1955) commonly observed in nature (Alves et al. 2001; Lada et al. 2003; Teixeira et al.2005), using the FLASH AMR code. The refinement criteria resolves the Jeans’ length ( λ J = p πc s /Gρ ) by at least eight cells, limited by the ambipolar diffusion time step. The ambipolar2iffusion time step is τ AD = 0 . x η AD = 0 . µ γ AD ρ i ρ n ∆ x . B , (1)where ∆ x is the cell length, η AD is the diffusivity, γ AD = 3 . × − cm s − is the ion-neutralcollisional coupling, B is the magnetic field strength, and ρ i = 2 . × − (cid:0) ρ n / − (cid:1) . g cm − (e.g. Duffin & Pudritz 2008).We construct an initial BE sphere with a radius of R = 6 . r = 6 . c s core / √ πGρ c (where thecritical radius is R c = 6 . r ) and add a 10% over–density on top with a 10%, m = 2 perturba-tion to break axisymmetry. Such perturbations model the buffeting that observed star–formingcores undergo from supersonic turbulent gas motions that prevail in the surrounding low densitymolecular gas (our perturbation accounts for a 0.04% change in gravitational binding energy, sothis is a rather quiescent core). Our BE spheres have a mass of 1 . M ⊙ and initial central densi-ties of ρ c = 6 . × − g cm − . The magnetic field is added to the core in the z direction – inagreement with recent observations (Kirk et al. 2009) – such that the ratio of the thermal to themagnetic energy in the gas, known as the plasma beta ( β = 2 c s /v = 46.01) is constant, where c s = 0 .
27 km s − is the core’s sound speed and the Alfvn speed is v A = B/ √ πρ ≃ . c s / Γ (thelatter relation for critical BE spheres). The mass-to-flux ratio Γ is a measure of the ratio of thegravitational to the magnetic energy in gravitating objects. If Γ <
1, the magnetic field is strongenough to prevent gravitational collapse and the core is called sub-critical, otherwise it is calledsupercritical. We set the mass-to-flux Γ = 2 πG / Σ /B = 3 .
5, where Σ is the column density,stemming from observations of magnetic fields in early cores which indicate Γ ≈ − β rot = 0 .
3. Do disks form?
By the time the maximum surface densities in the moderate rotation model set reach a valueof Σ c = 4 . × g cm − , 0.1 M ⊙ is contained within 156 AU (about 10% of the core’s mass iscontained within 0.003% of its initial volume). The mass of the protostar (taken to be all materialinside 1 AU or 10 − M ⊙ ) is only 6 % of the total mass in the disk (taken to be everything inside 10AU or 10 − M ⊙ ). At this point, the non–magnetic, perfectly coupled and partially coupled cloudsare 0.160, 0.194, and 0.175 million years old respectively (corresponding to 5.9, 7.2 and 6.5 free-falltimes, where t ff = 0 .
027 Myr is estimated using ρ c ). An extended dense region has taken shape,3he gas has been heated, and in the cases where magnetism is present, disk winds have started.We present in Figure 1, snapshots of the structure of the region at the center of the collapse,taken at this time for our three different cases. Despite having identical initial conditions, thenon–magnetized, moderate rotation case has formed two bars, driving off–axis rotational modeswhich themselves create an off–axis bar. The wobbling in the collapse occurs through accretionalong shocks and density perturbations in the gas. This stands in stark contrast to the obviousdisk–like structures and outflows present in both ambipolar diffusion (Figure 1(b) and ideal MHDcases (Figures 1(c) and (d) . Strong oscillations of the disk radius were noted in the simulations ofPrice & Bate (2007) which were also dominated by massive disks in their early phase. They alsofound that this behavior was suppressed by magnetic fields.Most importantly, the magnetic disks, by preventing strong bars from forming, are rotationallydominated whereas the uncoupled disk is not! This is illustrated in Figure 2(a), where we see thatthe ambipolar diffusion case has a larger rotationally dominated disk ( v φ > v infall ) than the perfectlycoupled case. We note also, by plotting v φ /v infall of the most evolved perfectly coupled state of thecollapse, that the rotational support is growing with time while hydrodynamic rotational supportbecomes increasingly unstable. In Figure 2(a) we also plot surface averaged plots of the ratio ofthe rotational velocity to the Keplerian velocity – the rotational velocity of a gas in a stable orbitaround an enclosed mass. We find peak ratios of about 0.1–0.5 in all cases of the moderate rotationmodel (in fact some pre–binned values exceed 1 in the ambipolar case).In the limit of low rotation, hydrodynamic spiral waves and bars in the purely hydrodynamiccase are confined to a region of about 1 AU in radius. This has little effect on the braking of thecloud (Figure 2(d) and rotationally dominated disks are allowed to form outside of 1 AU (Figure2(b). The hydrodynamic case is much more spherical in nature while the magnetized cases remainfairly flat, despite efficient magnetic braking. In the limit of low rotation, Figure 2(b) shows thatthe rotation in the hydrodynamic case is nearly Keplerian and the magnetic case is much moresub–Keplerian, whereas in the case of moderate rotation shown in Figure 2(a) the opposite is true.Thus in the extreme low rotation limit where spiral modes and bars are strongly suppressed, earlyhydrodynamic disks do form.Due to the bar, the distribution of specific angular momentum j ( R ) (where R is the cylindricalradius; Figure 2(d) of the hydrodynamic case is comparable to the cases where magnetic braking ispresent. In the low rotation model set, the bar is small and j ( R ) does not compare as well to themagnetized cases. The hydrodynamic and perfectly coupled cases of the moderate rotation modelset have similar maximal accretion rates of ˙ M accr ≈ × − M ⊙ yr − , while the ambipolar casereaches maximal accretion rates of only ˙ M accr ≈ × − M ⊙ yr − due to weaker magnetic brakingand a stable disk. In the low rotation model set, mass accretion rates are ˙ M accr ≈ × − M ⊙ yr − for the magnetized cases and ˙ M accr ≈ × − M ⊙ yr − for the hydrodynamic case. We provide movies of Figures 1(a)-(c) online illustrating the density and magnetic field line structure on all scalesin the simulation.
4n the high rotation model set however, no cases form rotationally dominated disks as all threecases have formed some sort of bar on large scales (Figure 3). Clearly, it is through the suppressionof instabilities in the collapse that magnetic fields promote early disk formation. Disk formation insystems without magnetic fields typically must wait until the mass of the central star dominatesthat of the surroundings sufficiently to suppress these instabilities.
4. Early outflows and magnetic fields
Figures 1(b)-(d) are snapshots of the structure of outflows that are launched at these earlytimes. A major result, shown by comparing Figures 1(b) and (c), is that outflows occur even whenambipolar diffusion is active. Outflows in a non–ideal MHD collapse have been demonstrated in two-dimensional ambipolar diffusion without drift heating (Mellon & Li 2009) and in three-dimensionalusing an ohmic diffusion approximation (Machida et al. 2007). Our analytical understanding of theearly outflow mechanism is the magnetic tower (Lynden-Bell 2003). In this picture, the toroidalfield in the disk is generated through the rotating flow and confined there by the accretion rampressure. As the field winds up, it reaches a critical strength in which the toroidal pressure is strongenough to overcome the ram pressure, moving wrapped field lines away from the disk and takingperfectly coupled gas with them – starting with the initial release in Figure 1(c) and the moredeveloped magnetic bubble shown in Figure 1(d). This acts to torque down the disk, removingangular momentum and increasing accretion rates . In the ambipolar diffusion picture, these fieldlines can seep through the ram pressure as they are not so tightly coupled to the accreting gas.We have also discovered the existence of a centralized, high–speed component to the ambipolardiffusion case’s outflow ( v > . − in an region roughly 1 AU in size). This componentresembles a jet that is commonly observed in the evolved perfectly coupled collapse . We note thatthe base of this component to the outflow lies on a cushion of heat generated by the friction betweenions and neutrals as toroidal field lines seep through the accreting gas. This spike in drift heatingis typical in ambipolar diffusion mediated C-shocks, regions where ion–neutral friction is strong(Wardle 1991). This local heating effect produces strong gradients in temperature with T >
T > B z is identical for all low and moderaterotation, magnetized cases (Figure 2(f), the toroidal field B φ in the low rotation, magnetized casesis nearly 2 orders of magnitude smaller than in the moderate rotation, magnetized cases.In Figure 2(c) we plot the ratio of mass loss rates ˙ M wind to mass accretion rates ˙ M accr for A movie of the development of Figure 1(d) from Figure 1(c) is provided online. The velocity structure of the perfectly coupled outflow is seen in our online movie of the evolution of Figure 1(d) − and observed outflow rates of ˙ M wind / ˙ M accr ∼ . yr due to time stepconstraints. However, from this early footprint, we are led to believe that the outflow mechanismis working. In the low rotation, ideal MHD case, the outflow rates are significantly diminishedby nearly 2 orders of magnitude in comparison to mass accretion rates, less than what is typicallyobserved. This may have consequences on feedback effects in star forming clusters and the assemblyof massive stars (Wang et al. 2009); protostellar energy injection is dependent on the rotationalenergy of the star forming core.
5. Fragmentation and binaries
In supercritical cores, three-dimensional ideal MHD calculations have shown that strong per-turbations (on the order of 50%) are needed in order for cores with Γ ∼ c = 68 g cm − ). The hydrodynamic case quickly fragments into a binary with a wideseparation of roughly 1000 AU. The perfectly magnetized case forms a bar and does not fragment,stabilized through its magnetic pressure as well as through the effects of magnetic braking. Mean-while, the ambipolar diffusion case produces an intermediate result. Early on a ring is formed whichfragments into a bar–like object with a small companion at about 1000 AU.For sufficiently high rotation ambipolar diffusion allows binaries to form without the need ofstrong (e.g. &
10% in amplitude) perturbations. For moderate rotation rates however, fragmenta-tion can take place only if large perturbation amplitudes ( ≈ . Conclusions In the early stages of a purely hydrodynamic collapse of a moderately rotating core, materialjoins the protostar by accreting through a chaotic series of bars and spiral waves. Our resultsshow that in magnetized collapses however, the magnetic field suppresses these wave modes, and asmall, regular disk appears at the earliest times. The weakening of magnetic control by ambipolardiffusion is insufficient to guarantee the formation of binary stars in typical cores with moderaterotation. Our results suggest that modest turbulent amplitudes ( > ∝ r − (1 . − . much more quickly with disk radius than do protoplanetary disk models at latertimes, wherein Σ ∝ r − . or r − . As the collapse winds up the field early outflows appear – evenfor partially coupled disks – and feed angular momentum and mechanical energy back into the starforming neighborhood. The density and magnetic structure of ideal and partially coupled disks arequite similar – the main difference is in the strength of the wound up magnetic field (at least anorder of magnitude weaker in the ambipolar diffusions case).Finally, we note that evidence for early localized intense drift heating near the disk at theaccretion shock, which heats materials up to 1,000K and beyond in a localized region, may bepreserved in the observed composition of comets and meteorites. The rapid heating, and subsequentcooling of those crystalline Mg–rich silicate materials that passed through the accretion shock andthrough the region of high ambipolar heating is reminiscent of heating events that must haveoccurred for some of these materials seen in cometary grains (Wooden et al. 2007). These eventshave been attributed to shock heating by spiral waves out to 10 AU in disks (Desch et al. 2005).As seen in Figure 1(b), only a portion of the accreted disk material will pass through this localizedheated region. Subsequent radial turbulent mixing of this flash heated material with the bulk ofmaterial that passed through a more gentle heating environment could in principle contribute tothe wide mixture of thermal histories preserved in cometary grain materials.
7. Acknowledgements
We thank Robi Banerjee, Eve C. Ostriker, and James Wadsley for fruitful discussions. D.D. wassupported by McMaster University and R.E.P. by the Natural Sciences and Engineering ResearchCouncil of Canada. We are pleased to acknowledge the SHARCNET HPC Consortium for the useof its facilities at McMaster University. The software used in this work was in part developed by theDOE–supported ASC/Alliances Center for Astrophysical Thermonuclear Flashes at the Universityof Chicago. This research was supported in part by the National Science Foundation under grantPHY05–51164.
REFERENCES
Alves, J. F., Lada, C. J., & Lada, E. A. 2001, Nat, 409, 1597rce, H. G., Shepherd, D., Gueth, F., Lee, C.-F., Bachiller, R., Rosen, A., & Beuther, H. 2007, inProtostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 245–260Banerjee, R., & Pudritz, R. E. 2006, ApJ, 641, 949—. 2007, ApJ, 660, 479Banerjee, R., Pudritz, R. E., & Anderson, D. W. 2006, MNRAS, 373, 1091Basu, S., & Mouschovias, T. C. 1994, ApJ, 432, 720Bonnor, W. B. 1956, MNRAS, 116, 351Crutcher, R. M. 1999, ApJ, 520, 706Desch, S. J., Ciesla, F. J., Hood, L. L., & Nakamoto, T. 2005, in Astronomical Society of thePacific Conference Series, Vol. 341, Chondrites and the Protoplanetary Disk, ed. A. N.Krot, E. R. D. Scott, & B. Reipurth, 849–+Duffin, D. F., & Pudritz, R. E. 2008, MNRAS, 391, 1659Ebert, R. 1955, Zeitschrift fur Astrophysik, 37, 217Hartmann, L., & Calvet, N. 1995, AJ, 109, 1846Hennebelle, P., & Fromang, S. 2008, A&A, 477, 9Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25Hosking, J. G., & Whitworth, A. P. 2004, MNRAS, 347, 1001Kirk, J. M., Crutcher, R. M., & Ward-Thompson, D. 2009, ApJ, 701, 1044Lada, C. J., Bergin, E. A., Alves, J. F., & Huard, T. L. 2003, ApJ, 586, 286Larson, R. B. 2009, ArXiv e-printsLynden-Bell, D. 2003, MNRAS, 341, 1360Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2007, ApJ, 670, 1198—. 2009, ApJ, 699, L157Mellon, R. R., & Li, Z.-Y. 2008, ApJ, 681, 1356—. 2009, ApJ, 698, 922Price, D. J., & Bate, M. R. 2007, MNRAS, 377, 77Pudritz, R. E., Ouyed, R., Fendt, C., & Brandenburg, A. 2007, in Protostars and Planets V, ed.B. Reipurth, D. Jewitt, & K. Keil, 277–2948eixeira, P. S., Lada, C. J., & Alves, J. F. 2005, ApJ, 629, 276Tilley, D. A., & Pudritz, R. E. 2007, MNRAS, 382, 73Wang, P., Li, Z.-Y., Abel, T., & Nakamura, F. 2009, ArXiv e-printsWardle, M. 1991, MNRAS, 251, 119Whelan, E. T., Ray, T. P., Bacciotti, F., Natta, A., Testi, L., & Randich, S. 2005, Nature, 435, 652Wooden, D., Desch, S., Harker, D., Gail, H.-P., & Keller, L. 2007, in Protostars and Planets V, ed.B. Reipurth, D. Jewitt, & K. Keil, 815–833
This preprint was prepared with the AAS L A TEX macros v5.2.
Model Set Low Rotation Moderate Rotation High RotationΩ t ff = Ω p π/ Gρ . . . − ) 1 . × − . × − . × − Rotational beta ( β rot ) 0 . .
046 0 . M core = 1 . M ⊙ , core temperature of T core = 20 K ( c s core = 0 .
267 km s − ), external temperature of T ext = 200 K ( c s ext = 0 .
845 km s − ),Γ = 3 .
5, and β = 46 .
01 and is run in three cases: hydrodynamic, ideal MHD, and ambipolardiffusion. 9ig. 1.— Inner disk of the moderate rotation model set at central surface densities of Σ =4 . × g cm − for (a) hydrodynamic, (b) ambipolar diffusion, (c) perfectly coupled cases, and (d)the evolved state of the perfectly coupled case at Σ c = 8 . × g cm − . Black and gray surfacesare ρ = 10 − g cm − ( n = 2 . × cm − ) and ρ = 10 − g cm − ( n = 2 . × cm − ),respectively. Temperature contours at T = 85 K are shown in blue. Outflow velocities contours,shown in red, are 0 . − in (c) and (d) and 0 . − . − in b. Magnetic field lines areshown as yellow tubes. (Animations of this figure are available in the online journal.)10 -2 -1 v φ / v K ep l e r R [cm]a) evolved idealevolved hydro10 -2 -1 v φ / v i n f a ll M od . R o t a t i on R [AU]a) idealnon-idealhydro 10 -2 -1 v φ / v K ep l e r R [cm]b) evolved idealevolved hydro10 -2 -1 v φ / v i n f a ll Lo w R o t a t i on R [AU]b) idealnon-idealhydro10 -6 -4 -2 Lo w R o t a t i on R [cm]c) idealnon-idealevolved ideal10 -6 -4 -2 -2 -1 d M / d t w i nd / d M / d t a cc r M od . R o t a t i on R [AU]c) 10 -5 -4 -3 -2 -1 Lo w R o t a t i on M enc [Msol]d)10 j [ c m s - ] M od . R o t a t i on d) j ∝ M enc1.00 j ∝ M enc0.93 idealnon-idealhydro10 Lo w R o t a t i on R [cm]e) 10 -1 Σ [ g c m - ] M od . R o t a t i on R [AU]e) Σ ∝ R -1.08 Σ ∝ R -2.51 Σ ∝ R -1.11 Σ ∝ R -1.73 idealnon-idealhydro 10 -2 Lo w R o t a t i on R [cm]f) 10 -2 B z , B φ [ µ G ] M od . R o t a t i on R [AU]f) B φ ∝ R -1.57 B z ∝ R -1.29 B z ∝ R -2.27 B φ ∝ R -1.47 B z ∝ R -1.20 B z ∝ R -1.67 idealnon-ideal Fig. 2.— Azimuthally averaged plots (density weighted, see Banerjee & Pudritz 2006) vs. diskradius for ideal MHD, non–ideal MHD, and hydro cases. In all figures the model sets are shownat a common central surface density of Σ c = 4 . × g cm − (moderate or mod. rotation) andΣ c = 3 . × g cm −2