Sub-Alfvenic Non-Ideal MHD Turbulence Simulations with Ambipolar Diffusion: II. Comparison with Observation, Clump Properties, and Scaling to Physical Units
aa r X i v : . [ a s t r o - ph . S R ] J u l Draft version October 29, 2018
Preprint typeset using L A TEX style emulateapj v. 11/10/09
SUB-Alfv´enic NON-IDEAL MHD TURBULENCE SIMULATIONS WITH AMBIPOLAR DIFFUSION: II.COMPARISON WITH OBSERVATION, CLUMP PROPERTIES, AND SCALING TO PHYSICAL UNITS
Christopher F. McKee
Physics Department and Astronomy Department, University of California, Berkeley, CA 94720; and Laboratoire d’Etudes duRayonnement et de la Mati`ere en Astrophysique, LERMA-LRA, Ecole Normale Superieure, 24 rue Lhomond, 75005 Paris, France
Pak Shing Li
Astronomy Department, University of California, Berkeley, CA 94720 andRichard I. Klein
Astronomy Department, University of California, Berkeley, CA 94720; and Lawrence Livermore National Laboratory,P.O.Box 808, L-23, Livermore, CA 94550
Draft version October 29, 2018
ABSTRACTAmbipolar diffusion is important in redistributing magnetic flux and in damping Alfv´en waves inmolecular clouds. The importance of ambipolar diffusion on a length scale ℓ is governed by theambipolar diffusion Reynolds number, R AD = ℓ/ℓ AD , where ℓ AD is the characteristic length scale forambipolar diffusion. The logarithmic mean of the AD Reynolds number in a sample of 15 molecularclumps with measured magnetic fields (Crutcher 1999) is 17, comparable to the theoretically expectedvalue. We identify several regimes of ambipolar diffusion in a turbulent medium, depending on theratio of the flow time to collision times between ions and neutrals; the clumps observed by Crutcher(1999) are all in the standard regime of ambipolar diffusion, in which the neutrals and ions are coupledover a flow time. We have carried out two-fluid simulations of ambipolar diffusion in isothermal,turbulent boxes for a range of values of R AD . The mean Mach numbers were fixed at M = 3 and M A = 0 .
67; self-gravity was not included. We study the properties of overdensities–i.e., clumps–in the simulation and show that the slope of the higher-mass portion of the clump mass spectrumincreases as R AD decreases, which is qualitatively consistent with Padoan et al. (2007)’s finding thatthe mass spectrum in hydrodynamic turbulence is significantly steeper than in ideal MHD turbulence.For a value of R AD similar to the observed value, we find a slope that is consistent with that of thehigh-mass end of the Initial Mass Function for stars. However, the value we find for the spectralindex in our ideal MHD simulation differs from theirs, presumably because our simulations havedifferent initial conditions. This suggests that the mass spectrum of the clumps in the Padoan et al.(2007) turbulent fragmentation model for the IMF depends on the environment, which would conflictwith evidence for a universal IMF. In addition, we give a general discussion of how the results ofsimulations of magnetized, turbulent, isothermal boxes can be scaled to physical systems. Eachphysical process that is introduced into the simulation, such as ambipolar diffusion, introduces adimensionless parameter, such as R AD , which must be fixed for the simulation, thereby reducing thenumber of scaling parameters by one. We show that the importance of self-gravity is fixed in anysimulation of ambipolar diffusion; it is not possible to carry out a simulation in which self-gravity andambipolar diffusion are varied independently unless the ionization is a free parameter. We show thatour simulations apply to small regions in molecular clouds, generally with ℓ . . M . M ⊙ .A general discussion of the scaling relations for magnetized, isothermal, turbulent boxes, includingself-gravitating systems, is given in the Appendix. Subject headings:
Magnetic fields—MHD—ISM: magnetic fields—ISM: kinematics and dynamics—stars:formation INTRODUCTION
Giant molecular clouds, threaded by magnetic fields,are the birth places for new stars. Since the earlieststudies of star formation, it has been recognized thatthe magnetic flux in stars is many orders of magnitudeless than that in the interstellar material from whichthe stars originated. Mestel & Spitzer (1956) suggestedthat ambipolar diffusion (AD) could resolve this prob- [email protected]@[email protected] lem by allowing magnetic flux to be redistributed duringcollapse due to the differential motion between the ion-ized and neutral gas. With effective shielding of highenergy cosmic rays and radiation, the ionization fractionof gas inside high-density cloud cores can be ≤ − (e.g.Caselli et al. 1998; Bergin et al. 1999), which rendersAD efficient. Star formation theory based on the AD-regulated, quasi-static collapse of molecular clouds (e.g.Spitzer 1968; Nakano & Tademaru 1972; Mouschovias1976, 1977, 1979; Nakano & Nakamura 1978; Shu 1983;Lizano & Shu 1989; Fiedler & Mouschovias 1992, 1993) McKee, Li, & Kleinnaturally accounts for the enormous loss of magnetic fluxduring star formation.However, both observations (Zuckerman & Evans1974; Zuckerman & Palmer 1974) and theory(Arons & Max 1975) have long indicated that supersonicturbulent motions are important in molecular clouds,and this turbulence has a major effect on star formation(Mac Low & Klessen 2004; Ballesteros-Paredes et al.2007; McKee & Ostriker 2007). The kinetic energy of thesupersonic motions is observed to be comparable to themagnetic energy of the clouds, so that molecular cloudsare in approximate equipartition (e.g. Crutcher 1999;Heiles & Troland 2005; Troland & Crutcher 2008). Itshould be borne in mind that the amplitude of turbulentfluctuations decreases with decreasing scale; for exam-ple, Goodman et al. (1998) and Barranco & Goodman(1998) find that the NH lines within ∼ . x —Mac Low et al. 1995), which is prohibitive athigh resolution (e.g. Nakamura & Li 2008). Li, McKee,& Klein (2006; hereafter LMK) developed the Heavy-IonApproximation, which takes advantage of the negligibleion inertia in regions of very low ionization and can accel-erate simulations of ambipolar diffusion by large factors.In the Heavy-Ion Approximation, the mass-weighted ion-ization is increased by a factor R ∼ and the ion-neutral coupling coefficient is decreased by the same fac-tor, so that the momentum transfer between ions andneutrals is unaffected. Using a semi-implicit two-fluidscheme proposed by Mac Low & Smith (1997) (see alsoT´ o th 1995), LMK tested the Heavy-Ion Approximationwith several classical problems involving ambipolar dif-fusion and found speed-ups of order a factor 100.In the first astrophysical application of the Heavy-IonApproximation, Li et al (2008, hereafter LMKF) stud-ied the statistical properties of supersonically turbulentsystems with ambipolar diffusion. The properties of theturbulence were found to vary smoothly from the hydro-dynamic case to the ideal MHD case as the importance of ambipolar diffusion decreased. They found that thepower spectra for the neutral gas properties of a stronglymagnetized medium with strong ambipolar diffusion aresimilar to those for a weakly magnetized medium; in par-ticular, the power spectrum for the neutral velocity isclose to that for Burger’s turbulence.In this paper, we extend this work on turbulent sys-tems with ambipolar diffusion but without self-gravity.This paper has three main goals: First, we give ageneral discussion of the ambipolar diffusion Reynoldsnumber, R AD , that characterizes ambipolar diffusion(Myers & Khersonsky 1995; Zweibel 2002) ( § R AD for the molecu-lar regions studied by Crutcher (1999) ( §
2) and showthat they are consistent with the theoretically expectedones. Second, we use numerical simulations to determinethe properties of the clumps that appear in a turbulentmedium with ambipolar diffusion ( §§ R AD . Third, we ana-lyze the scaling properties of simulations with ambipolardiffusion and determine the range of physical parame-ters that characterize the simulations ( § THE AMBIPOLAR DIFFUSION REYNOLDS NUMBER
The effects of ambipolar diffusion on a length scale ℓ in a medium with a flow velocity v can be character-ized by the ambipolar diffusion (AD) Reynolds number, R AD ( ℓ ). This quantity appears to have been first intro-duced by Myers & Khersonsky (1995); they referred to itas the magnetic Reynolds number, although that term isnormally used to describe the effects of Ohmic resistiv-ity. The AD Reynolds number is motivated as follows(Zweibel & Brandenburg 1997; Zweibel 2002, LMK):Ions in a partially ionized plasma are subject to twoforces: the Lorentz force, ∼ B / πℓ B , where B rms isthe rms magnetic field strength and ℓ B ≡ | B rms / ∇ B rms | ;and the drag force, γ AD ρ i ρ n v AD , where γ AD is the ion-neutral coupling coefficient, ρ i and ρ n are the ion andneutral densities, respectively, and v AD is the drift veloc-ity between the neutrals and the ions. When the ioniza-tion is low enough that the ion inertia can be neglected,these forces balance and the drift velocity is v AD ( ℓ B ) ≃ B πγ AD ρ i ρ n ℓ B . (1)We define the ambipolar-diffusion time over a lengthscale ℓ as t AD ( ℓ ) ≡ ℓv AD ( ℓ ) = 4 πγ AD ρ i ρ n ℓ B . (2)Similarly, we can introduce the ambipolar-diffusionlength scale ℓ AD , which is the length for which the am-bipolar drift velocity is the same as the flow velocity—i.e., in the frame of the ions, the length scale over whichHD Turbulence Simulations with Ambipolar Diffusion 3the field varies in a steady flow: ℓ AD = B πγ AD ρ i ρ n v . (3)In terms of the neutral-ion collision time, t ni = 1 /γ AD ρ i ,the AD time scale and length scale are t AD = ℓ v t ni , (4) ℓ AD = v t ni v , (5)where v A = B rms / (4 πρ ) / is the Alfv´en velocity andwhere we have assumed that the ion mass density is neg-ligible, so that ρ n ≃ ρ . The effect of ambipolar diffusionon a flow over a length scale ℓ with a characteristic ve-locity v is determined by the AD Reynolds number, R AD ( ℓ ) ≡ ℓvv A t ni = t AD t f = ℓℓ AD = 4 πγ AD ¯ ρ i ¯ ρ n ℓvB , (6)where t f ≡ ℓ/v is the flow time across a length ℓ . Ob-serve that ambipolar diffusion increases in importance as R AD ( ℓ ) decreases; thus, it becomes more important atlow densities, low ionizations, low velocities, small dis-tances and high field strengths. As Myers & Khersonsky(1995) showed, the ratio of the size of the region, ℓ ,to the minimum wavelength of a propagating Alfv´enwave in which the inertia is provided by the neutrals , λ min = πv A t ni (Kulsrud & Pearce 1969), is directly pro-portional to R AD ( ℓ ): ℓλ min = R AD ( ℓ ) π M A , (7)where M A ≡ v/v A is the Alfv´en Mach number.We have defined the AD Reynolds number in terms ofthe mean densities and the rms field strength. In a su-personically turbulent medium, the densities are subjectto large fluctuations, and if the Alfv´en Mach number islarge also, the magnetic field has large fluctuations aswell. If one defines R AD in terms of the local densitiesand field strength, then one can devise several differentways of averaging so as to obtain an effective value of R AD for a turbulent medium; in particular, the lengthscale ℓ can be taken to be the size of the region or it canbe determined self-consistently as the size of the aver-age eddy or density fluctuation. The resulting values for R AD in a turbulent region of size ℓ range from slightlylarger than R AD ( ℓ ) to several times less (see Paper IIIfor further discussion). One should thus bear in mindthat R AD ( ℓ ) is a characteristic value for the ratio of theambipolar diffusion time to the flow time, and the actualvalue in a turbulent medium might differ from this by afactor of a few.Mouschovias (private communication) has emphasizedthat the AD Reynolds number is useful for turbulentmedia in which the velocity dispersion is determined bythe turbulence (the case we are considering here), butnot in systems in which the flow velocity is determinedby the AD process itself. For example, in quasi-static,AD-regulated star formation, the AD length scale, ℓ AD is proportional to the radius of the self-gravitating cloud,and R AD is of order unity. Similarly, R AD is not a useful parameter to characterize C-shocks (Draine 1980), sincethe structure of such shocks adjusts itself so that R AD ∼ Numerical Evaluation of R AD ( ℓ )Evaluation of the AD Reynolds number requires eval-uation of both the ion-neutral coupling coefficient, γ AD ,and of the mean ionization mass fraction, ¯ χ i ≡ ¯ ρ i / ¯ ρ . Ifthis mass fraction is small ( ¯ χ i ≪ ρ n ≃ ¯ ρ and R AD ( ℓ ) = 12 (cid:18) γ AD ¯ ρℓc s (cid:19) ¯ χ i M β = (cid:18) γ AD ¯ ρℓc s (cid:19) ¯ χ i M A2 M , (8)where c s is the isothermal sound speed, M ≡ v/c s isthe Mach number and β ≡ π ¯ ρc s /B is the plasma β parameter. We normalize our results to the case in whichthe ionization is dominated by HCO + . The ion-neutralcoupling coefficient is then γ AD = 1 . × − cm s − m n + m i = 3 . × cm s − g − (9)(Draine, Roberge, & Dalgarno 1983), provided the rela-tive velocity of the ions and neutrals is less than about 19km s − . Note that this value of the coupling coefficientdiffers from that adopted in LMKF due to our assump-tion that the ionization is dominated by HCO + . Moregenerally, we shall write γ AD = 3 . × γ AD ∗ cm s − g − , (10)where γ AD ∗ is a number that allows for ions other thanHCO + .Next, we consider the ionization. The processes thatdetermine the ionization in molecular clouds are com-plex, and in general the ionization is time-dependent.We adopt a characteristic value of the ionization basedon the assumption that the ionization is in a steady stateand is dominated by HCO + . In equilibrium, the meanionization fraction by number is¯ x e, eq = (cid:18) ζ CR α ¯ n H (cid:19) / , (11)where ¯ n H = ¯ ρ/µ H is the mean density of H nuclei, µ H = 2 . × − g is the mass per H nucleus, ζ CR ∼ × − s − is the cosmic ray ionization rate perH atom (see the discussion by Dalgarno 2006) and α is the relevant recombination rate (McKee & Ostriker2007). Equation (11) is consistent with the results ofPadoan et al. (2004) at late times and high densitiesfor α = 2 . × − cm s − , the value they adoptedfor the dissociative recombination rate of HCO + . Ifsmall PAHs dominate the ionization, then the disso-ciative recombination rate is about 10 times smaller(Wakelam & Herbst 2008) and the ionization severaltimes larger. Tassis & Mouschovias (2005, 2007), whoincluded the effects of charged grains, adopted a disso-ciative recombination rate α = 1 . × − cm s − ; theirresults for the ionization are approximately consistentwith equation (11) for densities n H . cm − .Inference of the ionization from observations gener-ally requires knowledge of the cosmic ray ionizationrate and the density as inputs to the chemical models McKee, Li, & Kleinused to interpret the observations (Williams et al. 1998;Padoan et al. 2004). The ionization can be characterizedby the parameter C i defined by x e ≡ C i (cid:18) ζ CR n H (cid:19) / . (12)In the equilibrium model above, C i = α − / , whichis 630 cm / s / for the fiducial case. Williams et al.(1998) found a median ionization x e = 4 . × − (notethat they normalized their results to H , whereas we arenormalizing with respect to H). Their adopted ioniza-tion rate ( ζ CR = 2 . × − s − ) and density ( n H =5 × cm − ) correspond to C i ≃ / s / .More recently, Padoan et al. (2004) have interpretedthese data with time dependent models and infer lowervalues of the ionization and therefore C i . They find x e ≃ . × − − . × − and attribute the highervalues to the effect of FUV photoionization; for theirassumed ionization rate ( ζ CR = 6 × − s − ) and den-sity (2 × cm − ), the implied value of the ionizationparameter is C i ∼
600 cm / s / , fortuitously close toour fiducial value. The difference between the values of C i inferred by Williams et al. (1998) and Padoan et al.(2004) is a reflection of the uncertainties that remain indetermining the ionization in molecular clouds.To evaluate the AD Reynolds number, we require theion mass fraction ¯ χ i , which is related to the ion numberfraction ¯ x i ≡ ¯ n i / ¯ n H by ¯ χ i = ¯ x i m i /µ H → . x i , wherethe numerical evaluation is for HCO + . We then have¯ χ i = m i C i µ H (cid:18) ζ CR ¯ n H (cid:19) / = 2 . × − χ ∗ i ¯ n / , ! , (13)where the numerical factor χ ∗ i = m i
29 amu (cid:18) C i
630 cm / s / (cid:19) (cid:18) ζ CR × − s − (cid:19) / (14)allows for deviations from the fiducial case and where¯ n H , ≡ ¯ n H / (10 cm − ). Under the assumption thatthe ion mass is indeed about 29 amu (i.e., the mass ofHCO + ), the results of Williams et al. (1998) correspondto χ ∗ i ≃
3, whereas those of (Padoan et al. 2004) corre-spond to χ ∗ i ≃
1. The ionization can also be expressedin the form ρ i = C ρ / (Shu 1983), with C = χ i ρ / =1 . × − χ ∗ i g / cm − / . Shu (1983) adopted a value C = 3 × − g / , corresponding to χ ∗ i ∼
3, in agree-ment with the estimate of Williams et al. (1998). Nu-merically, the AD Reynolds number is then R AD ( ℓ ) = 16 . γ AD ∗ χ ∗ i M β ¯ n / , ℓ pc T / ! , (15)where T ≡ T / (10 K) and ℓ pc ≡ ℓ/ (1 pc). Regimes of Ambipolar Diffusion
We can distinguish several regimes in ambipolar diffu-sion in a turbulent medium. For χ i ≡ ρ i /ρ ≪
1, we have ρ n ≃ ρ so that the neutral-ion collision time, t ni and thecorresponding ion-neutral collision time, t in are related by t in ≡ γ AD ρ n = χ i γ AD ρ i = χ i t ni . (16)Similarly, the ion-Alfv´en Mach number, M A i ≡ v (4 πρ i ) / /B rms , is related to the Alfv´en Mach number, M A , by M A i = χ i M A2 . It follows that R AD = M A2 (cid:18) t f t ni (cid:19) = M A i (cid:18) t f t in (cid:19) . (17)We now identify five different regimes for ambipolardiffusion. For simplicity we ignore possible differencesbetween the velocity dispersions of the neutrals and ions(to be discussed in Paper III), which could change thecoefficient in front of M A i by up to a factor 2 in theexpressions below.I. Ideal MHD ( t f /t ni → ∞ , corresponding to R AD →∞ for a given value of M A ): The ions and neutralsare perfectly coupled.II. Standard AD ( t f > t ni ≫ t in , corresponding to R AD > M A2 ): The neutrals and ions are coupledtogether over a flow time so that the AD is weak.For M A = O (1), linear Alfv´en waves can prop-agate, since the propagation condition for Alfv´enwaves of wavelength ℓ derived by Kulsrud & Pearce(1969) is equivalent to R AD ( ℓ ) > π M A (eq. 7).The wave damping is weak (i.e., Γ t f <
1, whereΓ = k v t ni is the damping rate for low-frequencywaves) for R AD ( ℓ ) > π .III. Strong AD ( t ni > t f > t in , corresponding to M A2 > R AD > M A i ): The neutrals are no longercoupled to the ions in a flow time, but the ions re-main coupled to the neutrals. For M A = O (1),Alfv´en waves cannot propagate since λ min /ℓ >π/ M A & t in > t f , corresponding to M A i > R AD ): The ions and neutrals are onlyweakly coupled and act almost independently. Thedamping rate for the high-frequency Alfv´en wavesthat propagate in the ions is Γ = 1 / (2 t in ), sothese waves are weakly damped in this regime:Γ t f = R AD ( ℓ ) / (2 M A i ) < ℓ . The Heavy-Ion Approximation is basedon the assumption that the ion inertia is negligibleand therefore does not apply to this regime (seebelow).V. Hydrodynamics ( t f /t in → χ i →
0, corre-sponding to R AD → M A ):The neutrals are not affected by the trace ions andact purely hydrodynamically. One can of courserecover the hydrodynamic limit by letting B → M A → ∞ ; in that case, R AD is uncon-strained. It should be noted that the boundarybetween the hydrodynamic regime and the weaklycoupled regime is a matter of choice; if one de-mands that the ions have at most a 1% effect onthe neutrals, for example, then R AD would have tobe smaller than if one demands that the effects belimited to 10%.HD Turbulence Simulations with Ambipolar Diffusion 5It should be borne in mind that in all regimes exceptthe last (where it is irrelevant), we have assumed thatthe ions are well-coupled to the magnetic field—i.e., theion gyrofrequency is much larger than the ion-neutralcollision frequency, Ω i t in ≫
1. Although we have definedthe AD regimes for arbitrary values of the Alfv´en Machnumber (provided B is large enough that Ω i t in ≫ M A ∼ O (1), as it generally is in molecular gas inthe interstellar medium. Observed Values of R AD Crutcher (1999) has summarized sensitive Zeemanmeasurements of magnetic field strengths together withother physical parameters, including the plasma β andthe Mach numbers, for 27 molecular clouds. Of these,12 have only an upper limit on the line-of-sight magneticfield. Table 1 lists the values of the parameters fromTables 1 and 2 in Crutcher (1999) that we use to com-pute the corresponding R AD using equation (15). Wetake the length scale ℓ to be the cloud diameter. Weuse Crutcher’s correction for projection effects on themagnetic field: Zeeman observations determine the line-of-sight component of the field, B los , and on averagethe value of B that enters the plasma- β parameter is3 B . We assume that the parameters describing theion-neutral coupling and the ionization ( γ AD ∗ and χ ∗ i )are unity. From Table 1, we see that clouds with mea-sured field strengths have R AD ranging from a few to ∼
70. Because the range of R AD is so large, we quotethe logarithmic average, defined as h R AD i log ≡ h log R AD i ; (18)the logarithmic mean and dispersion of the AD Reynoldsnumber in these clouds is h R AD i log = 17 ± . h R AD i log = 22; if we discard L889 as an outlierbecause of its unusually high Mach number ( M = 7 . M A , in Table 1. All the cloudshave R AD > M A2 , implying that these clouds are inthe standard AD regime ( § −
15 for clouds with measured field strengths.We conclude that the effects of ambipolar diffusion mustbe considered in studies of molecular clouds, at least inthose regions shielded from the interstellar radiation fieldso that χ ∗ i = O (1), in agreement with studies extendingback for many years (e.g., Mouschovias 1987). Predicted R AD and Implied Self-Gravity As we now show, it is possible to predict the ADReynolds number for self-gravitating clouds that havean ionization of the form given in equation (13). As acorollary, we show that the strength of self-gravity is nota free parameter in simulations of ambipolar diffusion ina turbulent medium.The importance of self-gravity in a cloud of radius R or in a simulation box of size ℓ = 2 R is determined by the virial parameter (Bertoldi & McKee 1992), α vir ≡ σ R GM = 5 σ ℓ GM , (19)where σ = 1 √ M c s (20)is the 1D velocity dispersion in the cloud. The virialparameter is thus proportional to the ratio of kinetic togravitational energy. We wish to treat both real clouds,which we approximate as effectively spherical, and tur-bulent boxes. Of course, real clouds are not spherical(Bertoldi & McKee 1992 give the generalization to ellip-tical clouds), but keeping track of these two cases pro-vides a gauge of the importance of geometric effects; fur-thermore, the spherical cloud model has long been in use(e.g., Solomon et al. 1987). Let the area and volume ofthe cloud or box be A ≡ c A ℓ , (21) V ≡ c V ℓ , (22)where c A = ( π/ ,
1) and c V = ( π/ ,
1) for a sphericalcloud and a box, respectively. The virial parameter thenbecomes α vir = 5 M c s c V G ¯ ρℓ . (23)Since ¯ ρ i ≡ ¯ χ i ¯ ρ , equation (8) for the AD Reynolds numbercan be rewritten as R AD ( ℓ ) = (cid:18) c V Gα vir (cid:19) / γ AD ¯ χ i ¯ ρ / M A2 , (24)which shows that the AD Reynolds number is determinedby the ionization, the Alfv´en Mach number and the virialparameter. Insofar as the ionization is a function of thedensity, R AD ( ℓ ) will also depend on density. However,in the case of greatest interest, in which ¯ χ i ∝ χ ∗ i ¯ ρ − / ,where χ ∗ i is a number that is unity in the fiducial case(eq. 14), the AD Reynolds number is fixed at R AD ( ℓ ) = 19 . γ AD ∗ χ ∗ i (cid:18) π c V α vir (cid:19) / M A2 . (25)Molecular cloud cores and clumps with measured mag-netic fields are typically self-gravitating, with α vir ∼ M A ∼ γ AD and χ i , such regions have R AD ( ℓ ) ∼
20. This predicted value is in good agree-ment with the observed values discussed in § R AD ( ℓ )since most of their mass is photoionized by UV radiation(McKee 1989), so that they have a higher ionization thanthe cores and clumps within them (e.g., in an envelopeof a GMC in which the ionization is dominated by C + ,the ionization parameter is χ ∗ i ∼ ).The importance of self-gravity in a magnetized mediumcan also be expressed in terms of the ratio of the mass tothe magnetic critical mass, M Φ , which is the minimummass that can undergo gravitational collapse. In terms McKee, Li, & Kleinof the magnetic flux, Φ ≡ Bc A ℓ , the magnetic criticalmass is M Φ = c Φ Φ G / , (26)where c Φ = 1 / π for a cold sheet (Nakano & Nakamura1978) and ≈ .
12 for a cloud with a flux-to-mass dis-tribution corresponding to a uniform field threading auniform spherical cloud (Mouschovias & Spitzer 1976;Tomisaka et al. 1988). For c Φ = 1 / π , the ratio of themass to the magnetic critical mass is µ Φ , ≡ M M Φ = (cid:18) πc V c A α vir (cid:19) / M A (27) → (cid:18) (cid:19) / M A α / (spherical cloud) , (28)which provides a simple relation between the Alfv´enMach number, M A and the two parameters describ-ing the importance of self gravity in a magnetized, tur-bulent cloud, α vir and µ Φ , . The ratio µ Φ , is some-times written as the ratio of the observed mass-to-flux ratio to the critical one, ( M/ Φ) obs / ( M/ Φ) crit (e.g.,Troland & Crutcher 2008). Using equation (25), we findthat the AD Reynolds number is given in terms of µ Φ , by R AD ( ℓ ) = 13 . (cid:18) c A c V (cid:19) γ AD ∗ χ ∗ i M A µ Φ , ; (29)the factor in parentheses is unity for a spherical cloud.Gravitationally bound clouds that are both magnetizedand turbulent have µ Φ , somewhat greater than unitysince the gravity has to overcome both the turbulent mo-tions and the magnetic field (McKee 1989). This expres-sion thus gives a similar result to that in equation (25)for M A ∼ µ φ, ≃ − R AD ( ℓ ) can be inverted to give thevalues of the virial parameter and the ratio of the mass tothe critical mass in terms of R AD ( ℓ ) and M A . In otherwords, a simulation of a turbulent box with ambipolardiffusion [which requires specification of R AD ( ℓ ) and M A ] necessarily implies the strength self-gravity wouldhave were it to be included: α vir = 203 (cid:20) γ AD ∗ χ ∗ i R AD ( ℓ ) (cid:21) M A4 , (30) µ Φ , = 0 . (cid:20) R AD ( ℓ ) γ AD ∗ χ ∗ i (cid:21) M A , (31)where we have set c A = c V = 1, as is appropriate for asimulation box. For α vir ≫
1, the neglect of self-gravityis self-consistent. Parameter choices that lead to valuesof α vir ≪ µ Φ , & α vir & χ ∗ i is spec-ified. By contrast, a simulation of a turbulent box with Of course, the relation between α vir and R AD ( ℓ ) also depends ideal MHD is scale free; the density can be chosen sothat self-gravity would be negligible if it were included.This freedom does not exist in simulations of ambipolardiffusion. SIMULATIONS
In this paper, we extend the LMKF study of super-sonic turbulence with ambipolar diffusion, focusing onthe physical properties of the clumps formed purely asthe result of turbulent fragmentation with no gravity.LMKF performed a series of 256 simulations in a peri-odic box using the code ZEUS-MPAD to investigate tur-bulence statistics in non-ideal MHD without self-gravity.Like LMKF, we drove the turbulence with a fixed driv-ing pattern over the wavenumber range 1 ≤ k ≤ k ≡ k phys ℓ / π = ℓ /λ is the normalized wavenumber)using the recipe described in Mac Low (1999). The driv-ing maintained the 3D Mach number at M = 3, which isonly mildly supersonic. The corresponding line-of-sightMach number—i.e., the 1D Mach number M / √ β parameter of 0.1, corre-sponding to an Alfv´en Mach number M A = 0 .
67; theturbulence is thus sub-Alfv´enic. During the simulations,the volume-averaged magnetic field changed by less than10%, and as a result the volume-averaged value of β re-mained within 10% of its initial value. As shown in Table1, this value of β is close to the median of the 15 cloudswith measured magnetic fields.The focus of our effort is to determine how the prop-erties of the clumps vary with R AD ( ℓ ), so for now wediscuss our results in dimensionless form; the physicalconditions corresponding to these simulations will be dis-cussed in § ℓ ≃ . R AD ( ℓ ) from 0.12, close to thehydrodynamic limit, to 1200, close to the ideal MHDlimit. The run with R AD ( ℓ ) = 12 has conditions similarto those in observed clouds; as we shall see below, if weassume that the simulated region satisfies the linewidth-size relation, its density would be ¯ n H ≃ cm − . In § ℓ / − ℓ /
20, so the ADlength scale ℓ AD = ℓ /R AD ( ℓ ) is in the inertial range forthis run. For all the other runs, the AD length scale isoutside the inertial range. The run with R AD ( ℓ ) = 1 . ℓ AD and has a lower valueof the AD Reynolds number than any of the clouds ob-served by Crutcher (1999), most of which are gravitation-ally bound. If the simulation satisfied the linewidth-sizerelation, it would have a density of ¯ n H ≃
400 cm − , corre-sponding to an unbound cloud. On the other hand, theruns with R AD ( ℓ ) = 120 , ℓ AD and have higher AD Reynolds numbers thanany of the clouds with measured magnetic fields in thatsample. The R AD ( ℓ ) = 1200 run could be applied to theouter parts of GMCs, where the ionization is dominatedby C + . The run with R AD ( ℓ ) = 0 .
12 represents the on M A , but this is fixed in simulations of turbulent boxes with M A .
1; by contrast, whereas observed clouds have definite valuesof χ ∗ i , it is not necessary to specify this quantity in the simulation—see § A.3.
HD Turbulence Simulations with Ambipolar Diffusion 7transition to the hydrodynamic limit, and is primarily oftheoretical rather than practical interest.All models were run for 3 t f , where t f ≡ ℓ / M c s isthe flow time. In order to improve the statistics and theresolution, we re-ran the five models m3c2r-1 [ R AD ( ℓ ) =0 .
12] to m3c2r3 [ R AD ( ℓ ) = 1200; see Table 2] in LMKFwith the same initial conditions but using a 512 grid.All the results reported in this paper are the result ofsimulations on such a grid. The total computing time forall the models was ∼ χ i = R χ i and a corresponding ion-neutral cou-pling coefficient ˜ γ AD ∝ γ AD / R , with R ∼ ; herethe tilde denotes quantities measured in code units (see § A.3). The key to the Heavy Ion Approximation is thateven though each of these parameters differs from theactual value by a factor of 10 , the ion-neutral couplingis governed by the product of the parameters and hasthe correct physical value. According to the discussionin § R AD ( ℓ ) in equilibrium.Our second principal approximation was in our treat-ment of the ionization. Simulations can be carried outwith various assumptions about the ionization, includ-ing ion conservation, ionization equilibrium and time-dependent ionization. Following LMKF, we assumedthat the number of ions is conserved, so that the value of¯ χ i for the entire box is constant. The density is initiallyuniform, so that the initial ionization mass fraction, ¯ χ i ,is the same everywhere; we took it to be 10 − . LMKFdemonstrated that the results were the same as in thecase of ionization equilibrium (basically because the timefor a neutral to exchange momentum with an ion is smallcompared to the ionization time scale). More generally,the ionization is time-dependent. The ratio of the flowtime, t f ≡ ℓ / M c s , to the characteristic ionization time, t ion , eq = x e, eq /ζ CR (see eq. 11), is large: t f t ion , eq = ( αζ CR ¯ n H ) / (cid:18) ℓ M c s (cid:19) , (32)= 1 . × ( α ∗ ζ ∗ CR ) / (cid:20) R AD ( ℓ ) γ AD ∗ χ ∗ i (cid:21) M A2 , (33)where α ∗ ≡ α/ (2 . × − cm s − ), ζ ∗ CR ≡ ζ CR / (3 × − s − ), and we have used equation (15). It followsthat the molecular gas is typically very close to ionizationequilibrium (although it is not necessarily close to chemi-cal equilibrium). In simulations, the relevant comparisonis between the flow time across a cell, t f / N g , where N g is the number of grid cells in the length of the box, andthe ionization time. For our runs, which typically have N g = 512, we have t f / ( N g t ion , eq ) ≃ R AD ( ℓ ) for fidu-cial values of the parameters. Ionization equilibrium isthus a good approximation for all the cases we considerexcept R AD ( ℓ ) = 0 . models with time-dependent ionization fordifferent values of R AD . We find that the properties of the clumps in these runs are within a few percent of thosein the corresponding 256 runs with ion conservation,with the exception of the ion density. In fact, the meanion density in the entire box in the time-dependent caseis less than that in the ion conservation case by up to afactor ∼
2. As a result, the value of the AD Reynoldsnumber is reduced by a corresponding factor, as shown inTable 2. For large R AD ( ℓ ), the gas is close to ionizationequilibrium, so that ¯ ρ i ∝ ρ / . With this relation for theion density, the mean ion density, and hence R AD ( ℓ ),are reduced by only a small amount compared to thecase of ion conservation for the low Mach number weare considering if the density PDF is a lognormal witha width similar to that found by Padoan & Nordlund(2002). For small R AD ( ℓ ) the deviations from ioniza-tion equilibrium are larger, and correspondingly the dif-ference between the time-dependent and ion conservationresults are larger as well. In this paper, however, we areexploring the effects of changing R AD ( ℓ ) by orders ofmagnitude, so changes of . PHYSICAL PROPERTIES OF CLUMPS
The formation of high-density clumps is a naturaloutcome in simulations of highly supersonic turbulence,whether a magnetic field is included or not. Fur-thermore, high-resolution turbulence simulations (e.g.Li et al. 2004; Padoan et al. 2007) produce a mass spec-trum of clumps that qualitatively resembles the stel-lar initial mass function (IMF), with a peak at lowmass and a power-law tail at high masses. Recentobservations of molecular cores (e.g. Tachihara et al.2002; Onishi et al. 2002; Alves, Lombardi, & Lada 2007)suggest a similarity between the stellar IMF and thecore mass function. (We follow the terminology ofWilliams, Blitz, & McKee 2000 and use the term “core”to refer to the subset of clumps that are gravitation-ally bound and will form a star or small multiple stellarsystem.) Padoan & Nordlund (2002) and Padoan et al.(2007) have proposed a turbulent fragmentation theoryfor the IMF that relates the index of the velocity powerspectrum to the slope of the higher-mass end of theclump mass spectrum. LMKF showed that ambipolardiffusion changes the velocity power index, and we con-firm that conclusion in Paper III. If the turbulent frag-mentation theory is correct, we would expect a changein the slope of the higher-mass end of the clump massdistribution between the ideal MHD and the AD tur-bulence simulations as well. (It should be noted thatthe Hennebelle & Chabrier 2008 theory leads to a muchsmaller predicted difference in the slope of the IMF inthese two cases.)We use a CLUMPFIND algorithm, based on the algo-rithm developed by Williams, De Geus, & Blitz (1994),to determine the clumps in our simulations. We de-fine “clumps” as connected regions with a density largerthan the mean density of the turbulent box and willuse the term “ClMF” for “clump mass function,” reserv-ing “CMF” for “core mass function.” This distinctionis appropriate for our simulations since they do not in-clude self gravity. The density contours are separated by δρ = 0 . ρ , which Padoan et al. (2007) found to workwell in distinguishing distinct clumps. In order to inferthe effects of AD on the ClMF, we require the clumps McKee, Li, & Kleinto be resolved. As mentioned in LMK, ZEUS-MPADneeds at least 3 to 6 zones to accurately distinguish theeffects of AD from those of numerical diffusion. There-fore, we require clumps to have at least 6 zones in themean radius, unless otherwise specified; this requirementis validated in the resolution study of ClMF in § r c ≡ (3 V c / π ) / , where V c is thevolume of the clump is determined by summing the vol-umes of each cell in the clump that has a density abovethreshold; thus, for a porous clump, r c is less than theprojected radius of the clump (see § δρ , we found that thenumber of clumps with mean radius larger than 6 cellsdoes not change when the separation of the density con-tours is smaller than 4%, thereby justifying our choiceof δρ . Before constructing the ClMF, we verify that theclumps defined in our simulations satisfy the heavy-ionapproximation. The Heavy-Ion Approximation for Clumps
The condition for the validity of the heavy-ion approxi-mation is R AD ≫ M , where the ion Alfv´en Mach num-ber, M A i , is smaller than the total Alfv´en Mach number, M A , by a factor ( ρ i /ρ ) / ≪ R AD , c , we use the 3Ddensity-weighted velocity dispersion of the neutral gas, √ σ n , inside a clump as the flow velocity and the meandiameter of the clump, d c = 2 r c , as the length scale. Theion Alfv´en Mach number of the clump, M Ai , c , is takento be the rms value of M Ai of all the cells in the clump.We can re-write the definition of R AD in equation (6) forclumps as R AD , c ( D c ) ≡ πγ AD ρ i ρ n D c √ σ n B = γ AD ρ n D c σ n √ σ i M , c ≡ C HIA M , c . (34)In Figure 1, we plot M i,c versus R AD , c for modelsm3c2r-1, m3c2r1, and m3c2r3 at t = 3 t f ; the results formodels m3c2r0 and m3c2r2 lie between the nearby mod-els. The data points all have R AD , c ≫ M i,c , even formodel m3c2r-1, which has the smallest value of R AD ( ℓ ).We have verified that this is true at other times as well.LMKF found that the Heavy Ion Approximation wasvalid for a turbulent box provided R AD ( ℓ vi ) / M & ℓ vi is the length scale for ion-velocity variations,which is generally significantly smaller than the size ofthe box. We do not know how ℓ vi in the clumps com-pares with the clump diameters. If we assume thatthe two length scales are comparable, then the require-ment for the validity of the Heavy Ion Approximation is C HIA &
30. This is well satisfied for all the clumps ex-cept those in model m3c2r-1, which has R AD ( ℓ ) = 0 . C HIA ≃
10, and the Heavy Ion Approximationis at best marginally satisfied. We have not observed anyproblems associated with this, however.Two interesting features of the results are worth not-ing. First, almost all the clumps have smaller valuesof M A i than the box as a whole; this is expected be- cause of the linewidth-size relation. The few data pointswith slightly higher values of M A i are due to large sta-tistical fluctuation in the ion density in a few clumps.Second, we note that the distribution of the data pointsis roughly parallel to the power law R AD ∝ M A i (thestraight line). This is because the factor C HIA dependson two quantities, the column density, ρ n D c , and the ve-locity dispersion ratio, σ n /σ i , each of which is almostindependent of M A i . Clump Mass Function (ClMF)
Resolution: The Sonic Length and the Inertial Range
In studying the properties of the clumps that arisein boxes with supersonic turbulence, two length scalesare important: the sonic length, ℓ s , and the minimumscale for the inertial range, ℓ in , min , which correspondsto the wavenumber k in , max = ℓ /ℓ in , min . The soniclength, which is defined by the condition that the rmsturbulent velocity in a box of size ℓ s equal the soundspeed, gives a characteristic scale for density fluctu-ations in a supersonically turbulent medium (Padoan1995; V´azquez-Semadeni et al. 2003). The sonic lengthshould be well resolved in numerical simulations sinceit is important to resolve these density fluctuations andthe turbulent motions that produce them. The resolu-tion condition is ∆ x ≪ ℓ s , where ∆ x is the size of agrid cell; equivalently, in terms of the sonic wavenumber k s ≡ ℓ /ℓ s , we have ℓ / ∆ x ≡ N g ≫ k s . We assumethat the turbulence in the box exhibits a linewidth-sizerelation of the form M = 3 / σ nt c s = (cid:18) ℓ d ℓ s (cid:19) q , (35)where ℓ d is the effective minimum driving scale; the corre-sponding wavenumber is k d ≡ ℓ /ℓ d . In our simulations, k d = 2, and we find that the average Mach number inboxes of size ℓ d is indeed very nearly equal to that forthe entire box, M ≃
3. We also find q ≃ for R AD ( ℓ )in the range 0.12-12; for R AD ( ℓ ) = 120 , q ≃ . The sonic length in a simulation is then ℓ s = ℓ k d M /q . (36)Correspondingly, we have k s ≡ ℓ ℓ s = k d M /q ≃ − , (37)for q = and q = , respectively. This satisfies theresolution condition N g = 512 ≫ k s for R AD ( ℓ ) ≤ R AD ( ℓ ) = 120 , q = (the observed value–Heyer & Brunt 2004), thesonic length is related to the linewidth-size parameter Note that Krumholz & McKee (2005) defined the sonic lengthwith respect to the 1D turbulent velocity, σ nt = c s ( ℓ /ℓ s, ) q , andadopted q = ; the two versions of the sonic length are related by ℓ s, = 3 / (2 q ) ℓ s , corresponding to ℓ s, = 3 ℓ s for q = . HD Turbulence Simulations with Ambipolar Diffusion 9 σ pc (eq. 53) by σ c s ℓ s , (38)or ℓ s = 2 c s σ pc = 0 . T σ ∗ pc2 ! pc . (39)We define the inertial range of the turbulence as therange of wavenumbers over which the power spectrumis a power law in k . In our simulations, this extendsover the range k in , max > k >
3, where k in , max ≃ simulations and ≃
10 for the 256 simula-tions reported in LMKF. For k > k in , max , numerical dis-sipation becomes increasingly important. Another wayof expressing this condition is that with ZEUS, numer-ical dissipation becomes important at about 1/10th theminimum wavenumber, k in , max ≃ . × ( N g / k s < k in , max , and this is satisfied for the 512 simula-tions with R AD ( ℓ ) ≤
12. Determining whether this con-dition is a general requirement for accurate simulationsof supersonic turbulence is beyond the scope of this pa-per. We note that this condition becomes increasinglydifficult to satisfy as the Mach number increases.Figure 2 shows the clump mass distribution for the caseof R AD ( ℓ ) = 1200 (close to ideal MHD), at resolutionsof 512 and 256 . We can make an approximate relationbetween the clump masses and wavenumbers by associ-ating a wavenumber k c ≡ ℓ /D c , where D c is the clumpdiameter, to each clump. The corresponding clump massis approximately M c M = 4 π (cid:18) ¯ ρ c ¯ ρ (cid:19) k c ) , (40)where ¯ ρ c is the average clump density. For the high-resolution run, the mean density of the clumps withinthe inertial range ( k c <
20) is ¯ ρ c = 2 . ρ (i.e., the meandensity is 2.6 times the minimum clump density). Thehigher-mass part of the ClMF appears to be a power law(this is justified in § M c ∼ − .
3, correspondingto k c ≃
30 = 1 . k in , max . In fact, the clumps with sucha mass have D c ∼ −
20 cells. This is similar to boththe maximum wavenumber in the inertial range and tothe sonic wavenumber, which are also shown in Figure 2,to within a factor of 2. In order to determine whethereither of these parameters is associated with the changein slope, we also plot the clump mass spectrum for thecorresponding 256 run, for which k in , max = 10 (ver-tical dashed line) is reduced by a factor 2 whereas k s is unchanged. The results are clear: The break in theclump mass spectrum in the low-resolution run occurs athalf the wavenumber as in the high-resolution one. Fur-thermore, there is no discernable effect associated withthe sonic wavenumber, although it would be desirable totest this conjecture with both higher resolution simula-tions and for higher Mach numbers than M = 3, thevalue in the present simulations. It therefore appearsthat the dominant effect in determining the deviation ofthe ClMFfrom a power law is the numerical dissipationthat sets in for wavenumbers k > k in , max . The results of our simulations can therefore address only the higher-mass portion of the ClMF, with a minimum diameter of12 cells. Figure 3 shows the 3D spatial distribution ofclumps, identified by CLUMPFIND with minimum di-ameter of 12 cells, from a snapshot of model m3c2r1. Implications for the Turbulent Fragmentation Modelfor the IMF
As remarked above, the similarity between the coremass function and the stellar IMF suggests that theIMF may be defined during the formation of cores in-side molecular clouds. In the turbulent fragmentationmodel of Padoan & Nordlund (2002) and Padoan et al.(2007), the mass distribution of cores (i.e., gravitation-ally unstable clumps) has the form d N core d ln m ∝ (cid:20) (cid:18) m + σ x √ σ x (cid:19)(cid:21) m − Γ , (41)where σ x is the dispersion of the density PDF. Thepower-law index of the core mass function at high masses,Γ, is related to the index, n v , of the velocity power spec-trum, P ( k ) dk ∝ k − n v dk , byΓ = 34 − n v (42)for the strong-field, ideal MHD case ( B ≥ B cr ), andΓ = 35 − n v (43)for the weak-field, ideal MHD case ( B < B cr ), which in-cludes the hydrodynamic case (Padoan et al. 2007); herethe critical magnetic field B cr is defined by the conditionthat the postshock gas pressure be comparable to thepostshock magnetic pressure. Hennebelle & Chabrier(2008) have introduced an improved theory for the IMF,but we cannot comment on the differences between theirresults and those of Padoan et al. (2007) since our sim-ulations do not include self-gravity. As noted above, theHennebelle & Chabrier theory predicts a much smallerdifference in the slope of the IMF between the magneticand non-magnetic cases than does the theory of Padoanet al.The simulations of Padoan et al. (2007) do not includeself-gravity. Based on the discussion at the beginning of § M = 10), their box size of ℓ = 6 pc is in goodagreement with the linewidth size relation in equation(53). However, they chose a density ¯ n H = 2 × cm − ,which results in a virial parameter α vir ≃ . R AD ( ℓ ),as discussed in § M A we haveadopted ( M A = 0 .
67) and for the fiducial values ofthe ion-neutral coupling coefficient γ AD and the ioniza-tion parameter χ i , the virial parameter of the box is α vir = 41 /R AD ( ℓ ) (eq. 30). This is unphysically lowfor R AD ( ℓ ) &
6: in nature, large values of the AD0 McKee, Li, & KleinReynolds number are accompanied either by large val-ues of M A (which is unlikely according to the resultsof Crutcher 1999) or by larger ionizations than impliedby χ ∗ i = 1. The core Bonnor-Ebert mass—that is, theBonnor-Ebert mass based on the turbulent pressure inthe ambient medium (eq. A50)—in our simulations is M BE , core = √ M M BE = 0 . α / T σ ∗ pc2 M ⊙ → . R AD ( ℓ ) M ⊙ , (44)where the last expression is for fiducial values of the pa-rameters. Of the five AD models in Table 2, only modelm3c2r1, with R AD ( ℓ ) = 12 (comparable to the observedvalues), yields a physically plausible Bonnor-Ebert mass.We therefore do not attempt to put our clump mass func-tion in physical units here (physical units are discussedin § n v couldbe different, varying from the hydrodynamic relation tothe MHD one as R AD ( ℓ ) goes from 0 to ∞ . Second,as shown in LMKF, the value of n v also depends on R AD ( ℓ ). In the pure hydrodynamic case, n v = 2 for su-personic turbulence, whereas in the MHD case the valueof n v is not precisely known and could depend on theplasma β .In Figure 4, we show the ClMFs of models m3c2r-1,m3c2r1, and m3i. LMKF showed that the density cor-relation between data sets at different times approacheszero in a time slightly less than t f . Therefore, in orderto build up the statistics, we use three data sets in eachmodel run, at t ≃ t f , t f and 3 t f . Adding all the clumpstogether to form a single data set, we use reduced χ fit-ting to determine the higher-mass slope of the core massfunction, Γ fit . We face two problems in determining Γ fit :First, we do not know the the range of masses to includein the “higher-mass” data, and second, we do not wantour answer to depend on the size of the bins used in bin-ning the data. We begin by dividing the data into 20logarithmically spaced mass bins. To address the firstproblem, we carry out fits beginning with only the threehighest-mass bins, and then steadily increase the num-ber of bins used in the fitting until the peak of ClMF isreached. Initially, the value of χ drops as the numberof bins increases, since more data are contributing to thedetermination of the slope. However, when the numberof bins is large enough that the ClMF begins to deviatefrom a power law, the reduced χ will increase. To ad-dress the second problem, we increase the number of binsfrom 20 to 40 in increments of 5 and adopt the value ofΓ fit with the smallest reduced χ from all five sets of fit-ting. Usually, the slopes corresponding to the minimumreduced χ from different total bin numbers are close toeach other. The resulting slopes are listed in Table 3.In view of the noise fluctuations in the higher-massrange of the ClMF, we have performed a two-sampleKolmogorov-Smirnov (K-S) test to determine whetherthis part of the ClMF can be fit with a power law. Themass range extends from the highest mass bin to the breakpoint determined by the χ fitting procedure de-scribed above. The null hypothesis is that the higher-mass end of the ClMF from the simulation has the samedistribution as a power law. Our results show that theK-S test on all five AD models and the ideal MHD modelfails to reject the null hypothesis at the 5% confidencelevel. The p -values of all the K-S tests with different bin-ning are between 0.49 and 0.97. We conclude that thehigher-mass portion of the ClMF is statistically consis-tent with a power law. For the ideal mhd case (modelm3i) and for R AD ( ℓ ) = 12 (model m3c2r1), the powerlaw extends over the entire inertial range. However, inthe limit of low R AD ( ℓ ) (model m3c2r-1), the power lawextends only over the upper half of the inertial range;higher resolution and/or more samples are needed to de-termine if the inertial range is consistent with a powerlaw in this case.With these 512 models, the clump statistics are ade-quate to demonstrate that the higher-mass slopes dependon R AD . If turbulent fragmentation is correct, this is nosurprise because LMKF found that the spectral indexesof the velocity power spectra also depend on R AD . Herewe draw on the results of Paper III, which gives more ac-curate values of the spectral index for the velocity of theneutrals, n vn ( k ), than LMKF (see Table 3). The trendof spectral index changing from an Iroshnikov-Kraichnan(Iroshnikov 1963; Kraichnan 1965) to Burgers spectrum(Burgers 1974) as one goes from large to small R AD isstill clear, as reported in LMKF.In the limit of ideal MHD (model m3i), the higher-mass slope is Γ fit = 1 . ± .
09 (see Figure 3), whichagrees quite well with the prediction Γ = 1 .
18 from equa-tion (42) with spectral index n v = 1 . ± .
05. Notethat Padoan et al. (2007) get somewhat different results( n v = 1 . . M = 10 and β =1, whereas we have M = 3 and β = 0 .
1. In their hydro-dynamic simulations, Ballesteros-Paredes et al. (2006)found that the shape of the ClMF depends on the Machnumber of the turbulence, consistent with our result. Ifthe shape of the ClMF is significantly affected by the flowconditions, then the Padoan et al. (2007) model wouldimply that the IMF depends on the environment, sinceregions of star formation do not all have similar physi-cal conditions. As Ballesteros-Paredes et al. (2006) pointout, this could be problematic in view of observationalsupport for an IMF that is approximately universal.For the model m3c2r1, which has R AD = 12, compara-ble to the observed value ( § fit = 1 . ± .
10, which is consistent with the Salpetervalue. As noted above, in the limit of low R AD (modelm3c2r-1), we are unable to fit the data with a powerlaw that extends over the entire inertial range; The slopefor the high-mass portion of the range for which a fit ispossible is Γ fit = 2 . ± .
14, which continues the trendthat the slope increases as R AD ( ℓ ) decreases. Since thisslope applies to only part of the inertial range, however,we are unable to check the validity of equation (43),which relates the slope of the ClMF to the velocity powerspectrum in the weak field case. Consistent with the re-sults of Padoan et al. (2007), this slope is significantlygreater than the Salpeter value of the higher-mass slope,Γ = 1 . R AD modelsis greater than in the high- R AD models, except at thehigher-mass end: The total number of clumps with D c >
12 cells from the three time snapshots in thequasi-hydrodynamic model m3c2r-1 is 2093 (Table 3),whereas it is 1033 in model m3i. The total mass ofclumps in model m3c2r-1 is ∼ . M , whereas it is ∼ . M in model m3i. On average, the mass perclump in the quasi-hydrodynamic model m3c2r-1 issmaller than that in the ideal MHD model, which is alsoconsistent with prior simulations.The turbulent fragmentation model for the IMF pre-dicts that the core mass function (CMF) (i.e., the massfunction of gravitationally bound clumps) is the same asthe clump mass function (ClMF) at high masses, and itis based on the assumption that the IMF is proportionalto the core mass function (the latter is predicted in thework of Matzner & McKee 2000). Padoan et al. (2007)emphasize that the predicted ClMF for hydrodynamicturbulence is much steeper than for ideal MHD turbu-lence, and our work confirms this. Our work shows thatthere is a continuous variation in the higher-mass slopeof the ClMF due to the effects of ambipolar diffusion,such that the fraction of stars born at high mass shouldincrease with R AD . Mass-To-Flux Ratios
Ambipolar diffusion plays an important role in the corecollapse process when the clump mass is less than or com-parable to the magnetic critical mass (eq. 26). Observa-tionally, only a limited number of cores have measuredmass-to-flux ratios due to the difficulty in making pre-cise Zeeman measurements. Furthermore, observationsgive only the line-of-sight values for the magnetic fieldand column density, so the value of µ Φ ∝ M/ Φ for anyparticular core is necessarily uncertain. From a studyof 34 dark cloud cores, Troland & Crutcher (2008) foundan average value of µ Φ ,c = 1 . − . Resolution Study
In this section, we check the convergence of the mass-to-flux ratios of the clumps by comparing the results fromthe 256 and 512 runs for model m3c2r3. To carry outthe resolution study, we consider only clumps that havea mass at least equal to the minimum mass of clumpswith r c ≥
12 cells in the 512 model; for the 256 run,this corresponds to r c & µ Φ ,c /µ Φ , )/( M c /M ) / ver-sus M c /M for both the 256 and 512 runs in Figure 5. Curve fitting shows that the slope of the 256 data is0 . ± .
03 and the slope of 512 data is 0 . ± .
02. Themean values of ( µ Φ ,c /µ Φ , )/( M c /M ) / are 1 . ± . . ± .
03 for the 256 and 512 models, respectively.We conclude that the mass-to-flux ratios of clumps in the512 model are converged. Effect of R AD ( ℓ ) on the Mass-to-Flux Ratio As discussed in § R AD ( ℓ )on the mass-to-flux ratios, but at the expense of con-sidering models that would be unphysical were grav-ity to be included: Equation (31) implies µ Φ , =0 . R AD ( ℓ ) / ( γ AD ∗ χ ∗ i ) for our simulations, which is inthe observed range only for the R AD ( ℓ ) = 12 case. Whatis of interest then is how the normalized values of themass-to-flux ratio vary with R AD . For example, the ra-tio of the mass-to-flux ratio for an individual clump, µ Φ ,c ,to that for the entire box is µ Φ ,c µ Φ , = M c B c πR c, ⊥ · B ℓ M = (cid:18) B B c (cid:19) Σ c Σ ≃ Σ c Σ , (45)where R c, ⊥ is the radius of the clump normal to the fieldthreading the clump, B c , and Σ is the mean surfacedensity for the turbulent box. For the cases we consider,the mean field in the clump is close to the mean fieldof the whole box since the relatively small value of theAlfv´en Mach number, M A = 0 .
67, leads to a relativelyuniform field, as discussed in § / Σ c is justthe number of clumps along a flux tube. Furthermore,since the density of the clumps is typically a few times thethreshold density (see below eq. 40) and is thus approxi-mately constant, it follows that the mass-to-flux ratio inthe clumps is proportional to the cube root of the clumpmass: µ φ,,c ∝ Σ c ∝ ρ c R c ∝ R c ∝ M / c . (46)We now use our simulations to determine whether am-bipolar diffusion in a turbulent medium affects the mass-to-flux ratio in clumps, even in the absence of self-gravity.In order to ensure that the clumps we study are in thehigher-mass, power-law regime of the clump mass distri-bution so that numerical effects are minimal, we choosea minimum clump mass that is above the threshold forthe higher-mass regime in all cases. This minimum masscorresponds to a minimum clump radius of r c = 6 cells.To determine how ambipolar diffusion affectsthe mass-to-flux ratio, we compute the value of( µ Φ ,c /µ Φ , )/( M c /M ) / for all clumps in each modeland plot the results in Figure 6. The mean values of µ Φ ,c /µ Φ , for the three models are also tabulated inTable 3 and are shown as the horizontal lines in Figure6. In all the models, the values of µ Φ ,c for the clumps aresmaller than µ Φ , for the whole box due to fragmentationalong flux tubes. This effect has been observed in otherMHD turbulence simulations (V´azquez-Semadeni et al.2005a; Tilley & Pudritz 2007). The typical value ofthe normalized mass-to-flux ratio, µ Φ ,c /µ Φ , ∼ . µ Φ ,c scales as M / c .We observe from Figure 6 and Table 3 that h µ Φ ,c i /µ Φ , R AD model to the small and moderate R AD models.This table also shows that the mean density ofclumps, h ρ c i , in the three models increases systemati-cally as R AD decreases. The dispersion in the valuesof µ Φ ,c /µ Φ , / ( M c /M ) / in Figure 6 shows a significantvariation as R AD decreases. The dispersions of mass-to-flux ratio (not mass-to-flux divided by M / ) are givenin Table 3. The dispersion for R AD ( ℓ ) = 12 is almosttwice that for R AD ( ℓ ) = 1200. The larger dispersion of µ Φ ,c and higher density of clumps at R AD ( ℓ ) = 12 thanat high R AD ( ℓ ) suggest that material can more easilycross magnetic field lines as R AD decreases. A furtherdecrease in R AD ( ℓ ) to 0.12 results in a higher density,but increased fragmentation of the clumps reduces thedispersion somewhat. We conclude that, even in the ab-sence of self-gravity, ambipolar diffusion has an effect onthe mass-to-flux ratios of clumps. Other Physical Properties of Clumps
In this section, we summarize a number of other phys-ical properties of the clumps as functions of R AD , c in Figure 7 by comparing the two models m3c2r-1[ R AD ( ℓ ) = 0 .
12, strong AD] and m3c2r3 [ R AD ( ℓ ) =1200, strong ion-neutral coupling], which represent thetwo extremes of ion-neutral coupling among our simula-tions. Figure 7 gives side-by-side plots of the normal-ized clump radii, r c /ℓ , the ion and neutral densities, ρ i,c /ρ and ρ n,c /ρ , magnetic energy density, U B,c /U B, ,clump mass, M c /M , and ionization mass fraction, χ i ,for clumps for the two models.Figure 7a shows that the normalized radii of the clumpsin the strong-coupling model (m3c2r3) are, on average,larger than those for the strong AD model (m3c2r-1).This is a result of more fragmentation in the strongAD case. The largest radius in m3c2r3 is about dou-ble that in m3c2r-1. Since we require clumps to havea radius larger than 6 cells, there is a sharp trunca-tion in the size distributions at r c /ℓ = 6 /
512 = 0 . R AD , c because of the strong coupling betweenions and neutrals. This is seen in all other properties aswell. Figures 7b and 7c show the normalized mean ionand neutral densities of the clumps. The sharp bottomedge in Figures 7a and 7c is the result of the densitythreshold ρ c ≥ ρ we chose in defining the clumps. Thevariations in ion density are much smaller in the strongcoupling case than in the strong AD case. This is alsoreflected in the ionization mass fraction in Figure 7f. Theionization mass fraction of clumps in the strong couplingmodel is about constant, but that of clumps in the strongAD model varies by almost 3 orders of magnitude. InFigure 7d, the magnetic field is barely perturbed by theturbulence in model m3c2r-1 because of weak coupling;the magnetic field energy density in the clumps, U B,c , isvery nearly the same as that for the whole box. Althoughthe magnetic field is perturbed more in model m3c2r3,most clumps have U B,c within 50% of that in the box.Figure 7e shows that the largest clumps in m3c2r3 aremore massive than the largest ones in m3c2r-1. This larger mass is due to a larger size, since the densities inthe two models are about the same, and can be under-stood as the result of magnetic suppression of fragmenta-tion, as discussed in § r c = 6 cells. Fromthis figure, we see that the global physical properties ofclumps scale smoothly from r c = 6 cells to the largestclump. PHYSICAL UNITS FOR SIMULATIONS OF TURBULENTBOXES WITH AMBIPOLAR DIFFUSION
The results of our simulations have been reported indimensionless form. How can they be converted tophysical values? A simulation of an isothermal, mag-netized, turbulent box is characterized by three dimen-sional parameters—the size of the box, ℓ , the mean den-sity in the box, ¯ n H , and the isothermal sound speed, c s —and two dimensionless ones—the 3D sonic Machnumber, M = 3 / σ nt /c s and the plasma- β parame-ter, β ≡ π ¯ ρc s /B (Ostriker, Gammie, & Stone 1999;Padoan & Nordlund 1999). Here ¯ n H is the mean densityof hydrogen nuclei, σ nt is the 1D nonthermal velocitydispersion and B rms ≡ h B i / is the rms magnetic field.In the absence of other physical processes, all these pa-rameters can be selected arbitrarily, although the valueof c s ∝ T / is tightly constrained for molecular clouds,which generally have temperatures in the range 10 −
20 K.Inclusion of a new physical process, such as ambipolardiffusion, introduces a new dimensional constant, in thiscase the ion-neutral coupling parameter, γ AD . Corre-spondingly, a new dimensionless parameter (in this case, R AD ) can be formed and the number of independent di-mensional parameters is reduced by one. For a givensound speed, there is thus one independent dimensionalparameter, such as the density, in simulations of ambipo-lar diffusion; such simulations are therefore scale free.Treatments of ambipolar diffusion require specification ofthe ionization, which in principle can introduce anotherdimensionless parameter that in turn would determinethe scale. However, as discussed in § A.3, the Heavy IonApproximation eliminates this constraint.Adoption of a linewidth-size relation, as is observed inmolecular clouds (Larson 1981), also reduces the numberof independent dimensional parameters by one. Hence,if an isothermal system satisfies a linewidth-size relationand is subject to ambipolar diffusion, then its velocityscale is set by the isothermal sound speed, c s ∝ T / ,and its size and mean density are determined by dimen-sionless parameters. In this case, a given simulation ap-plies to only one set of parameters describing the box.This is discussed further in the Appendix, which givesexplicit expressions for properties of turbulent boxes inthe general case, when they satisfy a linewidth-size re-lation, and for self-gravitating boxes. Here we presentthe scaling for our simulations of turbulent boxes withambipolar diffusion. General Scaling Relations
To determine how simulations of a turbulent box withambipolar diffusion can be scaled to physical systems, weuse equation (15) to solve for the size of the simulationbox, ℓ . We find that it is determined by the remainingtwo dimensional parameters (¯ n H and T ) along with fiveHD Turbulence Simulations with Ambipolar Diffusion 13dimensionless parameters [ R AD ( ℓ ), γ AD ∗ , χ ∗ i , β and M ]: ℓ = 0 . (cid:20) R AD ( ℓ ) γ AD ∗ χ ∗ i (cid:21) MM A2 (cid:18) T ¯ n H , (cid:19) / pc , (47)The flow time across the box, the mass in the box, andthe column density are then t f = ℓ M c s = 1 . × (cid:20) R AD ( ℓ ) γ AD ∗ χ ∗ i (cid:21) M A2 ¯ n / , yr , (48) M = ¯ ρℓ = 1 . × − (cid:20) R AD ( ℓ ) γ AD ∗ χ ∗ i (cid:21) M T / M A6 ¯ n / , M ⊙ , (49) N H = ¯ n H ℓ = 9 . × (cid:20) R AD ( ℓ ) γ AD ∗ χ ∗ i (cid:21) MM A2 (¯ n H , T ) / cm − , (50)Note that these scalings are preserved by the Heavy-IonApproximation, in which the ion mass fraction ( ∝ χ ∗ i )is increased and the ion-neutral coupling coefficient ( ∝ γ AD ∗ ) is decreased by the same factor. The strength ofthe magnetic field does not depend on R AD ( ℓ ), B = (4 πρc s ) / MM A = 3 . n H , T ) / MM A µ G . (51)As discussed in § n − / have an implicit value of thevirial parameter, α vir , given by equation (30). Actualphysical systems have α vir &
1, since violations of thisinequality lead to gravitational motions that raise α vir up to order unity. Hence, this sets a lower limit on theproduct of the ionization parameter and the AD couplingparameter for a given value of R AD ( ℓ ), γ AD ∗ χ ∗ i & R AD ( ℓ )14 . M A2 . (52)The lower limit on the ionization corresponds to the caseof gravitationally bound clouds and clumps discussed in § . /c / V ≃ . Scaling with the Linewidth-Size Relation
Most molecular gas in the Galaxy is observed to obeya linewidth-size relation σ nt = σ pc R / , (53)where R pc is the radius of the region measured in pcand typically σ pc ≃ .
72 km s − (McKee & Ostriker2007). The linewidth-size relation is quite general: itapplies to within a factor ∼ R = ℓ / M is the 3D Mach number, we find M = 3 / σ nt c s = (cid:18) ℓ c s (cid:19) / σ pc (1 pc) / = 4 . σ ∗ pc ℓ , pc1 / T / , (54)where σ ∗ pc ≡ σ pc .
72 km s − . (55)Falgarone & McKee (2010) have shown that this turbulence-dominated linewidth-size relation applies only when N H < N LWS = 1 . × σ ∗ pc2 cm − , (56)or, equivalently, when¯ n H < ¯ n LWS = 9 . × σ ∗ pc4 M T ! . (57)For larger values of the column density and density,the linewidth-size relation must take the effects of self-gravity into account. The resulting virialized linewidth-size relation has σ ∝ (Σ ℓ ) / and is equivalent to settingthe virial parameter equal to unity, α vir = 1 (Heyer et al.2009; see § A.2.1). The linewidth is greater than that inthe turbulence-dominated case due to the effects of selfgravity.When the turbulence-dominated linewidth-size rela-tion applies, so that N H and ¯ n H satisfy the inequalities inequations (56) and (57), then the size of the simulationbox is determined by equation (54): ℓ = 0 . M T σ ∗ pc2 ! pc . (58)With the aid of equation (47), one can then express thedensity in terms of the linewidth-size parameter, σ ∗ pc ,¯ n H = 9 . × σ ∗ pc4 M T ! (cid:20) R AD ( ℓ )14 . M A2 γ AD ∗ χ ∗ i (cid:21) cm − , (59)= ¯ n LWS (cid:20) R AD ( ℓ )14 . M A2 γ AD ∗ χ ∗ i (cid:21) , (60)where the factor in brackets is . & N H = N LWS (cid:20) R AD ( ℓ )14 . M A2 γ AD ∗ χ ∗ i (cid:21) (61)with the aid of equation (56). The mass correspondingto N LWS and ¯ n LWS —i.e., the maximum mass at whichthe turbulence-dominated linewidth size relation holds—is M LWS , which is given in equation (A41). The mass inthe simulation box is given in terms of M LWS by M = M LWS (cid:20) R AD ( ℓ )14 . M A2 γ AD ∗ χ ∗ i (cid:21) . (62)On the other hand, when the system being simulatedis self-gravitating, then α vir ∼ ℓ ∝ ¯ n − / and M ∝ ¯ n − / are smaller than in the turbulence-dominated case, whereas N H ∝ ¯ n / is larger. The gen-eral case is discussed in the Appendix, § A.2.3.
Physical Parameters for Simulations
We are now in a position to discuss the physical param-eters corresponding to our simulations. For simplicity,we shall assume that the temperature is T = 10 K andthat the linewidth-size parameter has its fiducial value,4 McKee, Li, & Klein σ ∗ pc = 1, corresponding to σ pc = 0 .
72 km s − . Themaximum column density for the turbulence-dominatedlinewidth-size relation is N LWS = 1 . × cm − , andsince the Mach number is M = 3, the correspondingmaximum density is n LWS = 1 . × cm − . We shallfocus on the four cases R AD ( ℓ ) = 1 . , , , R AD ( ℓ ) = 0 .
12 case was done to model thetransition to the hydrodynamic limit. Recall that theclouds in Crutcher (1999)’s observations have a logarith-mic mean value h R AD i log = 17, comparable to the valuein the R AD ( ℓ ) = 12 simulation. The R AD ( ℓ ) = 1 . R AD , andthe R AD ( ℓ ) = 120 simulation a somewhat larger value,than any of the clouds in that sample; however, it mustbe borne in mind that this sample by no means covers allthe types of molecular gas in the Galaxy. In particular,the R AD ( ℓ ) = 1200 run is relevant to the outer parts ofGMCs, where the ionization is dominated by C + .The virial parameter associated with a given value of R AD ( ℓ ) in our simulations is α vir = (cid:20) R AD ( ℓ )6 . γ AD ∗ χ ∗ i (cid:21) − (63)from equation (30). Since our simulations have M A =0 .
67, the constraint on the ionization set by the require-ment α vir & R AD ( ℓ ) . . γ AD ∗ χ ∗ i . (64)First consider the case in which R AD ( ℓ ) = 1 . γ AD ∗ = χ ∗ i = 1);the ionization constraint is then well satisfied. Thevirial parameter is α vir = 28 from equation (63), sothe self-gravity is negligible in the system being simu-lated. Equations (47), (49), (50) and (51) imply thatthe size of the system is ℓ = 0 . n − / , pc, the massis M = 0 . n − / , M ⊙ , the column density is N H =7 . × ¯ n / , cm − , and the magnetic field is B =14¯ n / , µ G. Much of the unbound molecular gas in theGalaxy satisfies the turbulence-dominated linewidth-sizerelation (Falgarone et al. 2009). If the simulated systemsatisfies this relation, then the density is ¯ n H = 370 cm − from equation (60), and correspondingly the size of thesimulation box is ℓ = 0 . M = 0 . M ⊙ ,the column is N H = 4 . × cm − , and the magneticfield is B = 8 . µ G.Next, consider the simulations with R AD ( ℓ ) =12 , γ AD ∗ χ ∗ i is as close aspossible to its fiducial value, then the inequality in equa-tion (64) becomes an equality, and equations (47), (49),and (50) imply ℓ = 1 . n − / , pc, M = 82¯ n − / , M ⊙ ,and N H = 4 . × ¯ n / , cm − . These conditions cor-respond to α vir = 1, so the systems are on the virializedlinewidth-size relation. As remarked above, the virializedlinewidth-size relation applies when the density and col-umn density are large, ¯ n H & ¯ n LWS = 1 . × cm − and N H & N LWS = 1 . × cm − . Correspondingly, the size of the system is ℓ . . M . M LWS = 25 M ⊙ , and the magnetic field is B & µ G.In sum, the systems we have simulated are relativelysmall, with ℓ . / ¯ n / , pc (eqs. 47 and 52) and M . / ¯ n / , M ⊙ (eq. 49). If the system being simu-lated lies on the linewidth-size relation, then its mass is M . M LWS = 25 M ⊙ . As shown in § A.2.3, this inequal-ity also holds if the system has a linewidth greater thanthat given by the linewidth-size relation. Similarly, equa-tion (47) shows that the size of the system decreases with σ ∗ pc ; as a result, if the system being simulated lies on orabove the linewidth-size relation, then its size is no largerthan the size corresponding to the turbulence-dominatedlinewidth-size relation, ℓ . . CONCLUSIONS
Ambipolar diffusion is a key process in molecularclouds since it redistributes magnetic flux and dampswaves. The importance of ambipolar diffusion in a tur-bulent medium on a length scale ℓ and velocity disperion v is governed by the AD Reynolds number R AD = ℓ/ℓ AD ,where ℓ AD = v t ni /v is the length scale over which themagnetic field must vary in order to have a drift velocity v between the neutrals and ions (Zweibel & Brandenburg1997; Zweibel 2002). [Note that R AD is useful in de-scribing ambipolar diffusion whenever the velocity fieldincludes a significant turbulent component; it is not use-ful for non-turbulent, AD-driven gravitational collapse(Mouschovias 1987) or C-shocks (Draine 1980), where R AD is of order unity.] We have carried out two-fluidsimulations of isothermal, turbulent boxes using the codeZEUS MPAD (described in LMK) at a resolution of 512 for AD Reynolds numbers ranging from R AD = 0 .
12 to R AD = 1200, plus a simulation with ideal MHD. Theresolution we have used is sufficient to resolve the soniclength within the inertial range, permitting accurate sim-ulations for our M = 3 calculations. The mean Machnumbers were fixed at M = 3 and M A = 0 .
67, corre-sponding to a plasma- β parameter β = 0 .
1. The pur-pose of our simulations was to determine how the prop-erties of the clumps formed in molecular clouds dependon R AD ( ℓ ). One of the simulations (with R AD ( ℓ ) = 12)was in the middle of the observed range of the observedvalues of the AD Reynolds number; two of the simula-tions (those with R AD ( ℓ ) − . R AD ( ℓ ); andthe remaining two simulations, with R AD ( ℓ ) = 0 .
12 and1200, were designed to show the transition to hydrody-namics and ideal MHD, respectively. In order to carryout these simulations, we used the Heavy Ion Approxi-mation with an ionized mass fraction of 10 − (LMK) torepresent physical systems with actual ionized mass frac-tions ∼ − . We validated our simulations with conver-gence studies at lower resolution. The power spectra inour simulations show that the inertial range of our sim-ulations extends down to a length scale ℓ /
20, which iscomparable to the sonic length; it is important to resolvethe sonic length in simulations of turbulent boxes sincethe density has significant fluctuations on larger scales.Our principal conclusions are:HD Turbulence Simulations with Ambipolar Diffusion 151. Values of the AD Reynolds number R AD in a sam-ple of 15 molecular clumps with measured mag-netic fields (Crutcher 1999) range from 3 to 73; thelogarithmic mean value is 17. Omitting one out-lier, the clumps with upper limits on the magneticfield have an average lower limit of h R AD i log > R AD ≃
20, in excellentagreement with observation.2. Several regimes of ambipolar diffusion can be iden-tified, depending on the ratio of the flow time, t f , tothe ion-neutral collision time, t in , and the neutral-ion collision time, t ni : (I) ideal MHD ( t f /t ni → ∞ ,corresponding to R AD → ∞ for a given valueof M A ); (II) standard ambipolar diffusion, with t f > t ni , so that the neutrals and ions are cou-pled together over a flow time; (III) strong AD( t ni > t f > t in ), so that the neutrals are notcoupled to the ions over a flow time, but the ionsare coupled to the neutrals; (IV) weakly coupled( t in > t f ), so that the ions and neutrals behave al-most independently over a flow time; and (V) hy-drodynamics ( t f /t in → χ i →
0, correspondingto R AD → α vir that is de-termined by R AD , M A and the parameters de-scribing the ion-neutral coupling and the ionization[ α vir ∝ ( γ AD ∗ χ ∗ i /R AD ) M A4 —eq. 30]. It is notpossible to carry out a simulation in which the ef-fects of self-gravity and ambipolar diffusion are var-ied independently unless the ionization is treated asa free parameter.4. Clump mass spectrum. Using Clumpfind(Williams, De Geus, & Blitz 1994), we found allthe clumps with densities exceeding the mean den-sity in the box. We find that the slope of the higher-mass portion of the resulting clump mass spec-trum increases as R AD decreases, which is quali-tatively consistent with Padoan et al’s (2007) find-ing that the mass spectrum in hydrodynamic tur-bulence is significantly steeper than in ideal MHDturbulence. The value of the slope that we find for R AD = 12, the case closest to the value observedin molecular clouds, is Γ fit = 1 . ± .
10, whichis consistent with the Salpeter value, Γ = 1 . R AD = 1200) hasa slope Γ fit = 1 . ± .
11, which is marginallyconsistent with the Salpeter value. We furtherconfirm Padoan et al’s (2007) relation betweenthe index of the power spectrum and the slopeof the clump mass spectrum in the limiting casesof ideal MHD and near hydrodynamics. How-ever, the value we find for the spectral index inour ideal MHD simulation differs from theirs, pre-sumably because our simulation has lower values of β and M . This suggests that the IMF in thePadoan & Nordlund (2002); Padoan et al. (2007)turbulent fragmentation model depends on the en-vironment, which could conflict with evidence foran IMF that is approximately universal (see alsoBallesteros-Paredes et al. 2006).5. Ambipolar diffusion affects the mass-to-flux ratioof clumps, even in the absence of self-gravity: Theaverage mass-to-flux ratio µ Φ ,c /µ Φ , at low R AD isslightly larger than at high R AD , and the dispersionin the values of µ Φ ,c /µ Φ , for individual clumps atmoderate R AD is almost twice that at high R AD .6. Scaling relations for simulations of isothermal tur-bulent boxes. A simulation of an isothermal, mag-netized, turbulent box is characterized by threedimensional parameters: the size of the box, ℓ ,the mean density in the box, ¯ ρ ∝ ¯ n H , and thesound speed, c s (e.g., Ostriker, Gammie, & Stone1999). A single simulation with ideal MHD ap-plies to an infinite range of values of each of thesedimensional parameters, provided that the dimen-sionless parameters describing the simulation (inthis case, the Mach numbers M and M A ) are thesame (e.g., Padoan & Nordlund 1999). Except inregions of high-mass star formation, molecular gasgenerally has a temperature T ∼ −
20 K, so that c s is nearly constant; as a result, there are onlytwo dimensional parameters that have a significantvariation, ℓ and ¯ ρ . Each physical process thatis introduced into the simulation, such as ambipo-lar diffusion, introduces a dimensionless parameter,such as R AD , which must be fixed for the simula-tion, thereby reducing the number of scaling pa-rameters by one. For simulations with ambipolardiffusion, the physical parameters describing thesystem being simulated are characterized by a sin-gle dimensional parameter (for constant c s ), whichwe took to be the mean density ( § χ i ∝ ¯ n − / , the simulation remains scalefree. However, if one further requires that the sim-ulation satisfy an observed linewidth-size relation,then the mean density is determined and there areno independent scaling parameters.7. Physical parameters associated with the simula-tions. Two of the simulations we carried out werefor the purpose of studying the transition to idealMHD [ R AD ( ℓ ) = 1200] and to hydrodynamics[ R AD ( ℓ ) = 0 . R AD ( ℓ ) = 1 . , , α vir ∼
1. This con-straint gives a lower bound on the ionization suchthat R AD ( ℓ ) . . γ AD ∗ χ ∗ i . As a result, two ofthe simulations [ R AD ( ℓ ) = 12 , γ AD ∗ χ ∗ i = 1.The R AD ( ℓ ) = 1 . γ AD ∗ χ ∗ i = 1), whereas the R AD ( ℓ ) = 12 ,
120 sim-ulations correspond to systems that mostly likelyare gravitationally bound (for γ AD ∗ χ ∗ i as close tounity as possible). Because of this constraint andbecause of the small value of the Mach number weadopted ( M = 3), our simulations apply to smallregions in molecular clouds, with ℓ . / ¯ n / , pcand M . / ¯ n / , M ⊙ . If the system being sim-ulated has a velocity dispersion on or above thelinewidth-size relation observed in the Galaxy, thenthe size of the region is ℓ . . . M ⊙ .8. A general discussion of scaling relations for self-gravitating systems is given in the Appendix.In applying the linewidth-size relation, we follow Falgarone & McKee (2010) in distinguishing theturbulence-dominated relation from the virializedone.We thank Charles Hansen, Patrick Hennebelle, MarkHeyer, Mark Krumholz, Enrique Vazquez-Semadeni,Ellen Zweibel, and particularly an anonymous refereeand Telemachos Mouschovias for helpful comments. Thisresearch has been supported by the NSF under grantsAST-0606831 and AST-0908553 and by NASA under anATFP grant, NNX09AK31G. CFM also acknowledgesthe support of the Groupement d’Int´erˆet Scientifique(GIS) “Physique des deux infinis (P2I)” at the comple-tion of this work. RIK received support for this work pro-vided by the US Department of Energy at Lawrence Liv-ermore National Laboratory under contract DE-AC52-07NA 27344. This research was also supported by thegrant of high performance computing resources from theNational Center of Supercomputing Application throughgrant TG-MCA00N020. APPENDIX
SCALING LAWS FOR ISOTHERMAL TURBULENT BOXES
In this Appendix, we give a general discussion of scaling laws for simulations of isothermal, turbulent gases in a box.Although we do not include the effects of self-gravity in the text, we do include it here, so as to make the discussionmore generally useful. We focus on molecular gas, since such gas is generally approximately isothermal. Particularscaling relations that have been derived previously are noted (Klessen et al. 2000; Ostriker, Gammie, & Stone 1999;Tilley & Pudritz 2004; V´azquez-Semadeni et al. 2005a, 2008).In the simplest case in which there is no gravity and the MHD is ideal, a simulation of an isothermal, magnetized,turbulent box is characterized by two dimensionless parameters (Padoan & Nordlund 1999), the 3D sonic Mach number, M = 3 / σ nt /c s , and the plasma- β parameter, β ≡ π ¯ ρc s /B . Here σ nt is the 1D nonthermal velocity dispersion, c s the isothermal sound speed, ¯ ρ the mean mass density and B rms ≡ h B i / the rms magnetic field. Equivalently, thetwo parameters can be chosen to be the sonic Mach number and the Alfven Mach number, M A ≡ / σ nt /v A , sincethe plasma- β parameter is related to M and M A by β = 2 (cid:18) M A M (cid:19) . (A1)In general, the Mach numbers, M and M A , and the plasma- β parameter are functions of time.The turbulent box is also characterized by three dimensional parameters: the size of the box, ℓ , the meandensity in the box, ¯ ρ = M /ℓ , where M is the mass in the box, and the isothermal sound speed, c s (e.g.,Ostriker, Gammie, & Stone 1999). In the absence of other physical processes, these parameters can be selected arbi-trarily. In other words, a given simulation corresponds to definite values of M and β , but it can be scaled to arbitraryvalues of ¯ ρ , ℓ , and c s . However, the introduction of a new physical process, such as self-gravity or ambipolar diffu-sion, introduces a new dimensional constant and a corresponding new dimensionless parameter, so that the number ofindependent dimensional parameters is reduced by one. The same reduction occurs if a relation between dimensionalparameters is assumed, such as a relation between the size of the box and the mean velocity dispersion (a linewidth-size relation). In many cases, the temperature is tightly constrained, so that in fact there are only two dimensionalparameters that can be chosen at will. Hence, if an isothermal system satisfies a linewidth-size relation and is eitherself-gravitating or subject to ambipolar diffusion, then its velocity scale c s ∝ T / is set by the assumed temperatureand its size and mean density are determined by dimensionless parameters; in this case a given simulation applies toonly a single set of parameters describing the box.We return to the simplest case in which there is neither self-gravity nor ambipolar diffusion. Interstellar densitiesare often given in terms of number densities; we use the the density of hydrogen nuclei, n H = ρ/µ H , where µ H is themass per hydrogen nucleus (= 2 . × − g for cosmic abundances). Numerically, we have for the mass, flow time,and column density of the box, M = ¯ ρℓ = 34 . n H , ℓ , pc3 M ⊙ , (A2) t f ≡ ℓ v rms = ℓ M c s = 5 . × ℓ , pc M T / ! yr , (A3)HD Turbulence Simulations with Ambipolar Diffusion 17 N H = ¯ n H ℓ = 3 . × ¯ n H , ℓ , pc cm − , (A4)where T ≡ T / (10 K) and c s = 0 . T / km s − for molecular gas with cosmic abundances. The column densitycorresponds to a surface densityΣ = 2 . × − N H , g cm − = 11 . N H , M ⊙ pc − , (A5)where N H , ≡ N H / (10 H cm − ). The visual extinction corresponding to this column is A V = N H , δ mag,where δ is the ratio of the extinction per unit mass to the Galactic value. The magnetic field is given by(Ostriker, Gammie, & Stone 1999; note that their β is half the normal value) B rms = (cid:18) π ¯ ρc s β (cid:19) / = 4 . (cid:18) ¯ n H , T β (cid:19) / µ G , (A6)= 3 .
23 (¯ n H , T ) / MM A µ G . (A7) Scaling Relations for MHD Simulations of Turbulent Boxes with Self Gravity
As discussed above, self-gravity introduces an additional dimensionless parameter into a simulation and thereforereduces the number of independent dimensional parameters by one. For the case in which the ionization scales as¯ χ i ∝ ¯ n − / , this reduction is the same as that due to the inclusion of ambipolar diffusion ( § c A = c V = 1—see eqs 21 and22). With two dimensional parameters specified—the strength of self gravity and the temperature—there is still onefree dimensional parameter; as a result, self-gravitating, magnetized turbulent boxes are scale free. Dimensionless parameters
There are several equivalent dimensionless parameters that can describe the effects of self-gravity. One is the ratioof the mass to the characteristic mass of a self-gravitating cloud, c s / ( G ¯ ρ ) / : µ ≡ M c s / ( G / ¯ ρ / ) . (A8)This parameter is related to the mass, length and sound speed by GM /ℓ c s = µ / . (A9)In terms of the free-fall velocity, v ff ≡ ( GM / R ) / , we have µ = ( v ff /c s ) .Another parameter describing the effects of self-gravity is the number of Jeans lengths in the box (e.g.Ostriker, Gammie, & Stone 1999). The typical Jeans length in the box is λ J = ( πc s /G ¯ ρ ) / , so the number of Jeanslengths in the box is n J ≡ ℓ λ J = µ / √ π . (A10)The Jeans mass for the box is usually defined as M J ≡ ¯ ρλ , so that the number of Jeans masses in the box is M /M J = ( ℓ /λ J ) = n . The corresponding value of µ is µ J = M J ( G ¯ ρ ) / /c s = π / ≃ . R is α vir = 5 σ RGM (A11)(Bertoldi & McKee 1992); here σ = c s + σ is the total 1D velocity dispersion. Self gravity is important for α vir ≃ α vir ≫
1. By contrast, µ and n J can have arbitrary values & µ and n J , the effects of bulk kinetic energy as well as thermal energy are included in α vir . Forgas in a box, we define α vir by replacing R by ℓ / α vir ≡ σ ℓ GM , (A12)which is the same as equation (19) in the text. There is a complication here, since α vir is defined with respect to thetotal velocity dispersion, σ , whereas the linewidth-size relation depends only on the non-thermal velocity dispersion, σ nt . For M ≫
1, there is no problem, since the two velocity dispersions are nearly the same. Relations involving the8 McKee, Li, & KleinMach number that do not depend on the linewidth-size relation can be extended to low Mach numbers by redefining M as [3(1 + σ /c s )] / ; otherwise, such relations are restricted to M ≫
1. Bearing this in mind, we note thatequation (A9) implies that α vir is related to the other two parameters by α vir = 56 M µ / ! = 56 π (cid:18) M n J (cid:19) . (A13)V´azquez-Semadeni et al. (2008) derived a similar expression for spherical clouds (their α = α vir and their n J is R/λ J ,which is half the value we use). The virial parameter is also related to the ratio of the Jeans length to the sonic length( § α vir = 56 c s G ¯ ρ M ℓ s (A14)for q = . It follows that the ratio of the Jeans length to the sonic length is λ J ℓ s = (cid:18) πα vir (cid:19) / M . (A15)The parameters that describe the effects of self-gravity determine the ratio of the flow time to the free-fall time,which is t ff = (cid:18) π G ¯ ρ (cid:19) / = 1 . × ¯ n − / , yr . (A16)Relative to the free-fall time, the flow time is t f t ff = 1 . µ / M = 3 . n J M = 1 . α / . (A17) Variants of the Jeans length and Jeans mass
We can define both large-scale and small-scale variants of the Jeans length and Jeans mass. On large scales, thedensity is close to the mean, ¯ ρ , but the velocity dispersion is σ ≡ ( σ + c s ) / = M c s / √
3. We therefore define the“turbulent” variants of the Jeans length and Jeans mass by replacing the sound speed c s with the velocity dispersion σ , λ J , turb ≡ (cid:18) πσ G ¯ ρ (cid:19) / = (cid:18) M√ (cid:19) λ J , (A18) M J , turb ≡ ¯ ρλ , turb = (cid:18) M√ (cid:19) M J . (A19)The corresponding dimensionless quantities are µ turb , ≡ M σ / ( G / ¯ ρ / ) = (cid:18) √ M (cid:19) µ , (A20) n J , turb ≡ ℓ λ J , turb = (cid:18) √ M (cid:19) n J . (A21)Both µ turb , and n J , turb are of order unity when the virial parameter is: α vir = 52 µ / , = 52 πn , turb . (A22)On small scales, however, the velocity dispersion is about equal to the sound speed, c s , whereas the density can varyover orders of magnitude. In star-forming cores, the typical pressure is the mean turbulent pressure ¯ ρσ (Padoan 1995;Krumholz & McKee 2005); for an isothermal gas, this corresponds to a density ρ core ≡ ( M / ρ . We now introduceanother variant of the Jeans length, the “core Jeans length,” λ J , core , in which the velocity dispersion and density arethose expected in star-forming cores, λ J , core = (cid:18) πc s Gρ core (cid:19) / = (cid:18) πc s G M ¯ ρ/ (cid:19) / = (cid:18) √ M (cid:19) λ J . (A23)For supersonic flows, the core Jeans length is indeed small compared to the turbulent Jeans length, λ J , core =(3 / M ) λ J , turb , since it measures the effect of high pressures on thermally supported gas, whereas λ J , turb measuresHD Turbulence Simulations with Ambipolar Diffusion 19the effect of the turbulence on gas at the average density. If self-gravity is important in the turbulent box ( α vir ∼ λ J , core ℓ s = (cid:18) πα vir (cid:19) / = 3 . α / . (A24)The “core Jeans mass” is smaller than the normal Jeans mass and much smaller than the turbulent Jeans mass, M J , core = ρ core λ , core = 13 M ¯ ρλ , core = (cid:18) √ M (cid:19) M J . (A25)We expect M J , core (or perhaps the somewhat smaller core Bonnor-Ebert mass) to be the typical mass of gravitationallybound cores in a turbulent cloud (Padoan 1995). Scaling in Terms of the Mean Density
As discussed at the outset, the introduction of an additional physical process, such as self-gravity, reduces the numberof independent dimensional parameters to two, which we take to be the density n H and the temperature T . Since thetemperature has little variation in molecular clouds, there is effectively only one independent dimensional parameter, n H . In terms of n H , T and the dimensionless parameters describing the self gravity, the size and mass of the box arethen given by ℓ , pc = 0 . µ / (cid:18) T ¯ n H , (cid:19) / = 0 . n J (cid:18) T ¯ n H , (cid:19) / = 0 . M α / (cid:18) T ¯ n H , (cid:19) / , (A26) M M ⊙ = 4 . µ T / ¯ n / , ! = 22 . n T / ¯ n / , ! = 3 . M T / α / ¯ n / , ! , (A27) N H , = 1 . µ / (¯ n H , T ) / = 2 . n J (¯ n H , T ) / = 1 . (cid:18) ¯ n H , T α vir (cid:19) / M . (A28)Note that the column density is directly proportional to the square root of the thermal pressure in the first two cases,and to the square root of the turbulent pressure, ¯ ρv ∝ M ¯ n H T , in the third case; this is to be expected, since thepressure in a self-gravitating system is proportional to G Σ . The flow time is given by equations (A16) and (A17),the magnetic field by equation (A7) and the ratio of the mass to the magnetic critical mass, µ Φ , , by equation (27).The scaling relation between the size and density in terms of n J has been given by Ostriker, Gammie, & Stone (1999)and V´azquez-Semadeni et al. (2005a); for the mass in terms of the size and n J by Tilley & Pudritz (2004), althoughthey have a different numerical coefficient than implied by the above relations; and by Klessen et al. (2000) for boththe mass and the size in terms of the density for the particular case they consider, which has n J = 4.The Jeans mass and the Bonnor-Ebert mass are M J = 22 . T / ¯ n / , ! M ⊙ , M BE = 4 . T / ¯ n / , ! M ⊙ , (A29)where for the Jeans mass we have assumed that ¯ n H is the mean density in the ambient medium, and for the Bonnor-Ebert mass we have assumed that ¯ n H is the density at the surface of the Bonnor-Ebert sphere. Note that M J and M BE can be expressed in terms of the surface density of the box as, for example, by M J = 3 . µ / T N H 22 ! M ⊙ , M BE = 0 . µ / T N H , ! M ⊙ . (A30)We have changed the normalization of the column density so as to yield values of the Jeans mass and Bonnor-Ebertmass comparable to observed values for N H , ≡ N H / (10 cm − ) ∼
1. The core values are smaller by a factor √ / M .In this case, it is convenient to express the results in terms of the virial parameter of the box, M J , core = 5 . T α / N H , ! M ⊙ , M BE , core = 1 . T α / N H , ! M ⊙ . (A31)Note that in both cases, the critical masses have no explicit dependence on the Mach number.The basic conclusion is that MHD simulations of of self-gravitating, turbulent boxes are scale free; even with thetemperature fixed, there is one free parameter, such as the box size or the density, that can be chosen arbitrarily.Similarly, MHD simulations with ambipolar diffusion are scale free, as we have seen in §
5. If the ionization scales as1 / ¯ n / , then R AD and α vir are directly related ( § Scaling with the Linewidth-Size Relation
Linewidth-Size Relations for Molecular Clouds and Turbulent Boxes
Molecular gas in the Galaxy exhibits a linewidth-size relation in which the velocity dispersion of the gas increasesas a power of the physical dimension of the region (Larson 1981). Heyer & Brunt (2004) have shown that this applieswithin GMCs as well as among different GMCs. The data are consistent with the relation for the 1D velocity dispersion σ nt = σ pc R q pc , (A32)with σ pc = 0 .
72 km s − and q = (Solomon et al. 1987; McKee & Ostriker 2007); for these parameters, the Machnumber is M = 6 . R pc /T ) / . Such a relation appears to be satisfied by most molecular gas in the Galaxy;for example, Heyer & Brunt (2004) find q = 0 . ± .
15 within individual molecular clouds in their sample, andFalgarone et al. (2009) find that the combined data from several different surveys shows a clear linewidth-size relation,although they do not give a fit. The parameters have different values in regions of high-mass star formation, however:in such regions, σ pc ∼ a few km s − , different regions of high-mass star formation do not have line widths that increaseas R / , and it is not known how the velocity dispersion within individual regions scales with size (Plume et al 1997).Nonetheless, equation (A32) appears to be satisfied in these regions to within an order of magnitude.Mouschovias & Psaltis (1995) and Heyer et al. (2009) have proposed variants of the linewidth-size relation for gasthat is gravitationally bound. Falgarone & McKee (2010) have reconciled these virialized linewidth-size relations with the classical turbulence-dominated linewidth-size relation . Recall that the surface density of the cloud is Σ = M / ( c A ℓ ), where c A = (1 , π/
4) for a box and a spherical cloud, respectively. The virial parameter (eq. A12) is then α vir ≡ σ Gc A Σ ℓ , (A33)where we have used the identity symbol to emphasize that this follows directly from the definitions of the quantitiesinvolved; there is no physics in this relation. Solving this relation for the velocity dispersion gives σ ≡ (cid:20) π (cid:18) c A π/ (cid:19) Gα vir Σ R (cid:21) / , (A34)where R ≡ ℓ / c A is unity for a spherical cloud. This relationis superficially like the linewidth-size relation in equation (A32), but it is quite different: First, the exponent in thelinewidth-size relation, q ≃ , follows from observation, whereas that in equation (A34) is by definition, and secondthe coefficient in relation (A34) depends on the column density. Heyer et al. (2009) inserted the physics into thisrelation by noting that gravitationally bound clouds have α vir ≃
1. For spherical clouds, they then found σ = h π G Σ R i / , (A35)= 0 . N H , R pc ) / km s − (A36)which we term the virialized linewidth-size relation. They obtained a sample of bound clouds by combining the COdata on the Solomon et al. (1987) molecular clouds with data from the higher resolution CO data on these clouds fromthe Galactic Ring Survey (REF). Over a range of surface densities 10 M ⊙ pc − . Σ . M ⊙ pc − , equation (A36)describes the data well, after a somewhat uncertain correction is made for the cloud masses. Mouschovias & Psaltis(1995) previously found an analogous relation for magnetized clouds with Bµ Φ , in place of N H (see eq. 27).Falgarone & McKee (2010) concluded that non-self gravitating interstellar gas obeys the turbulence-dominatedlinewidth-size relation given by equation (A32), whereas self-gravitating gas satisfies the virialized linewidth-size rela-tion given by equation (A35). There is a critical surface density that defines the boundary between the turbulent andvirialized cases: Equating the velocity dispersions for the two cases givesΣ LWS = 5 πG (cid:18) π/ c A (cid:19) σ ! , (A37)= 192 (cid:18) π/ c A (cid:19) σ ∗ pc2 M ⊙ pc , (A38)where σ ∗ pc = σ pc / (0 .
72 km s − ) is normalized to the standard Galactic value (see eq. 55) and c A = π/ N LWS = 1 . × (cid:18) π/ c A (cid:19) σ ∗ pc2 cm − . (A39)It should be noted that these values for Σ LWS and N LWS for spherical clouds are comparable to the mean valuesfor Galactic GMCs found by Solomon et al. (1987), and about twice the mean values found by Heyer et al. (2009).HD Turbulence Simulations with Ambipolar Diffusion 21Regions with Σ ≪ Σ LWS are dominated by interstellar turbulence, whereas those with Σ ≫ Σ LWS are dominated byself-gravity and are decoupled from the turbulent cascade in the interstellar medium. The high surface densities ofregions of high-mass star formation thus naturally lead to the high velocity dispersions observed there by, for example,Plume et al. (1997). The maximum density for a cloud satisfying the turbulent linewidth-size relation is¯ n LWS = (cid:18) c A c V (cid:19) N LWS ℓ = 1 . × (cid:18) π/ c V (cid:19) σ ∗ pc4 M T cm − . (A40)where we used equation (54) to eliminate ℓ and the relations M ∝ c A N H ℓ = c V ¯ n H ℓ to cover the different geometries.The maximum mass for a cloud satisfying the turbulent linewidth-size relation is M LWS = 0 . M T σ ∗ pc2 ! M ⊙ , (A41)which is independent of geometry. Turbulent Boxes with the Turbulence-Dominated Linewidth-Size Relation ( N H ≤ N LWS , c A = c V = 1) Here we determine the scaling relations for turbulent boxes that satisfy the turbulence-dominated linewidth-sizerelation σ nt = 0 . σ ∗ pc R / km s − . (A42)The Mach number M and box size ℓ are related by equation (54). As a result, the properties of the turbulent boxare given by: ℓ = 23 M c s ( σ / . M T σ ∗ pc2 ! pc , (A43) t f = 2 . × M T / σ ∗ pc2 ! yr , (A44) M = 3 . × − M ¯ n H , T σ ∗ pc6 ! M ⊙ , (A45) N H = 1 . × M ¯ n H , T σ ∗ pc2 ! cm − . (A46)For simulations that include self-gravity but have surface densities less than the critical one (Σ ≤ Σ LWS ), the columndensity is most simply expressed in terms of the linewidth-size parameter and the virial parameter using equations(A33) and (A43), N H = 54 σ ! Gµ H α vir = 1 . × σ ∗ pc2 α vir ! cm − , (A47)which corresponds to N H = N LWS /α vir for c A = 1 in equation (A39). The scaling for the density can be expressed interms of the linewidth-size parameter σ ∗ pc and a parameter describing the self gravity with the aid of equations (A43)and (A26) , ¯ n H , = 115 µ / σ ∗ pc4 M T ! = 361 n σ ∗ pc4 M T ! = 96 σ ∗ pc4 α vir M T ! . (A48)Comparison with equation (A40) for a box geometry ( c V = 1) shows that ¯ n H = n LWS /α vir , where α vir &
1. Similarly,equation (A27) implies M M ⊙ = 0 . µ / M T σ ∗ pc2 ! = 1 . n M T σ ∗ pc2 ! = 0 . M T α vir σ ∗ pc2 ! , (A49)so that M = M LWS /α vir , where M LWS is given in equation (A41). The turbulence-dominated linewidth-size relationdoes not apply for densities exceeding n LWS , corresponding to column densities N H > N LWS and masses
M > M
LWS ,and for that case the scaling is given by the relations in § A.2.3 below.When expressed in terms of the linewidth-size relation, the core values for the Jeans mass and Bonnor-Ebert mass(eq. A25) are independent of the Mach number, M J , core = 3 . α / T σ ∗ pc2 M ⊙ , M BE , core = 0 . α / T σ ∗ pc2 M ⊙ . (A50)2 McKee, Li, & KleinThe core Bonnor-Ebert mass is comparable to the typical mass of observed stars, particularly if allowance is made forthe fact that only a fraction of the core mass is incorporated into the final star (e.g. Matzner & McKee 2000). Theserelations can be expressed in terms of the sonic length instead of σ ∗ pc by using equation (39).In applying these relations to simulations with driven turbulence, it must be kept in mind that the driving generallyresults in deviations from the linewidth-size relation (A32) on the driving scale. Using the linewidth-size relation torelate simulations to actual systems is therefore best done for cases in which the driving is restricted to large scales;the simulations discussed in the text satisfy this constraint since they are driven over a narrow range of wavenumbersat the largest scale, 1 ≤ k ≤ k d , with the driving wavenumber k d = 2. (Here k is a dimensionless wavenumber thatis related to the physical wavenumber k phys by k ≡ k phys ℓ / π ; the minimum possible wavenumber is k = 1 and themaximum is N g /
2, where N g is the number of grid cells in each side of the box.) When the turbulence is driven, onemust distinguish between the linewidth-size relation applied to the entire box, and the linewidth-size relation insidethe box. In our simulations, the mean Mach number is approximately constant over the range 1 ≤ k ≤ k d = 2. Wehave chosen to use the full size of the box in relating our simulations to clouds: ℓ = 2 R , where R is the radius of thecloud or of a region inside the cloud. However, in determining properties inside the cloud, such as the sonic length, itis necessary to allow for the fact that the internal linewidth-size relation is normalized approximately to the drivingscale. As a result, for σ ∝ ℓ / , the sonic length is ℓ s ≃ ( ℓ /k d ) / M rather than ℓ / M (see eq. 36). General Scaling with the Linewidth-Size Relation
The virialized linewidth-size relation follows from assuming that the virial parameter is unity, so the scaling relationsfor this case are given by the results in § A.1.3 with α vir = 1. To cover both the turbulence-dominated and virializedcases, note that equation (A46) shows that N H ∝ ¯ n H for the turbulence-dominated case ( α vir & n H . ¯ n LWS , whereas equation (A28) shows that N H ∝ ¯ n / for α vir = 1, corresponding to ¯ n H & ¯ n LWS . As a result, wehave N H = N LWS min " ¯ n H n LWS , (cid:18) ¯ n H n LWS (cid:19) / (A51)for the turbulence-dominated and virialized cases, respectively, as can be verified by direct substitution using equations(A28), (A39) and (A40). Similarly one can show that M = M LWS min " ¯ n H ¯ n LWS , (cid:18) ¯ n LWS ¯ n H (cid:19) / (A52)with the aid of equations (A27), (A40) and (A41). Note that M LWS is the maximum possible mass for a cloud witha given velocity dispersion, σ ∝ M T / , and linewidth-size coefficient, σ ∗ pc (Falgarone & McKee 2010). Furthermore,for a simulation with a given Mach number, M LWS decreases as σ ∗ pc increases (eq. A41); as a result, M LWS is also themaximum mass of a cloud with a linewidth above the linewidth-size relation. The size of the simulation box is ℓ = 0 . M T σ ∗ pc2 ! min " , (cid:18) n LWS ¯ n H (cid:19) / pc (A53)based on equation (A26). It must be borne in mind that these equations are based on the mean linewidth-size relation;for a given size and/or surface density, the velocity dispersion can vary by a factor of a few. Thus the virializedlinewidth-size relation, which applies to regions with high column densities by definition, correspondingly applies toregions of high density but with sizes and masses that decrease as the column density increases. Code Units
Numerical codes are generally written in dimensionless form, with masses, lengths and times written in terms ofcode units, ˜ M = M/M code , ˜ ℓ = ℓ/ℓ code , and ˜ t = t/t code . The code units can be adjusted to fit the problem beingsimulated. The properties of the box in code units, ˜ M = M /M code and ˜ ℓ = ℓ /ℓ code , can be selected arbitrarilyprior to the simulation (e.g., ˜ M = 8 and ˜ ℓ = 2), as can the normalized sound speed, ˜ c s = c s t code /ℓ code . If thereare N g grid cells in each side of the simulation box, the grid size is ∆˜ ℓ = ˜ ℓ / N g . For stationary gas, the time step is∆˜ t = C ∆˜ ℓ/ ˜ c s , where C is the Courant number.The code unit for length is given by ℓ code = ℓ / ˜ ℓ , where ℓ is given by equation (A43) if the typical Galacticlinewidth-size relation is adopted, by equation (A26) for a self-gravitating gas, and by equation (47) for a gas undergoingambipolar diffusion. The corresponding code unit for time is given by t code = ˜ c s ℓ code /c s . The code unit for mass isgiven by M code = M / ˜ M , where M is given by equation (A27) for a self-gravitating gas, by equation (A49) fora self-gravitating gas that obeys the linewidth-size relation, and by equation (49) for a gas undergoing ambipolardiffusion.The gravitational constant in the code is˜ G = GM code t ℓ = (cid:18) GM ℓ c s (cid:19) ˜ ℓ ˜ c s ˜ M = µ / ˜ ℓ ˜ c s ˜ M , (A54)HD Turbulence Simulations with Ambipolar Diffusion 23from equation (A9). Including the Heavy Ion Approximation ( ˜¯ χ i = R ¯ χ i and ˜ γ AD ∝ γ AD / R ), the ambipolar diffusionconstant in the code is ˜ γ AD = γ AD M code t code R ℓ = (cid:20) R AD ( ℓ ) R ¯ χ i M β (cid:21) ˜ c s ˜ ℓ ˜ M , (A55)where the second step follows from equation (8). [Keep in mind that ˜ ℓ , ˜ c s , and ˜ M are arbitrary; in the simulationsdescribed in the text, we have taken ˜ ℓ = 2, ˜ c s = 0 . ρ = ˜ M / ˜ ℓ = 1, so that ˜ γ AD = 0 . R AD ( ℓ ) / R ¯ χ i M β .] Solong as the Heavy Ion Approximation is valid, the outcome of a simulation is independent of the value of R since ˜ γ AD always enters in combination with ˜¯ χ i = R ¯ χ i .We can now address the issue of scaling in AD simulations when the physical ionization is specified by, for example,the value of χ ∗ i . In carrying out a simulation, the ionization in code units, ˜¯ χ i = R ¯ χ i , must be specified; hence, thephysical ionization χ ∗ i ∝ ¯ χ i ∝ R − . As noted above, the results of a simulation are independent of R so long as theHeavy Ion Approximation is valid. Thus, a single simulation provides the results for a family of problems with differentdegrees of ionization but the same values of R AD ( ℓ ) and c s ; as a result, we can use a single simulation to treat thephysically plausible range of ionizations for a given value of R AD ( ℓ ), as discussed in § REFERENCESAdams, F. C. & Shu, F. H. 2007, ApJ, 671, 497Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17Arons, J., & Max, C. E. 1975, ApJ, 196, L77Ballesteros-Paredes, J., Gazol, A., Kim, J., Klessen, R. S.,Jappsen, A.-K., & Tejero, E. 2006, ApJ, 637, 384Ballesteros-Paredes, J., Klessen, R. S., Mac Low, M.-M., &Vazquez-Semadeni, E. 2007, Protostars and Planets V, 63Barranco, J. A. & Goodman, A. A. 1998, ApJ, 504, 207Bergin, E. A., Plume, R., Williams, J. P., & Myers, P. C. 1999,ApJ, 512, 724Bertoldi, F. & McKee, C. F. 1992, ApJ, 395, 140Burgers, J. M. 1974, The Nonlinear Diffusion Equation(Dordrecht: Reidel)Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998,ApJ, 499, 234Crutcher, R.M. 1999, ApJ, 520, 706Dalgarno, A. 2006, PNAS, 103, 411Draine, B. T. 1980, ApJ, 241, 1021Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, 1983, ApJ,264, 485Duffin, D. F. & Pudritz, R. E. 2008, submitted to MNRAS,astro-ph 0810.0299Falgarone, E., Pety, J., & Hily-Blant, P. 2009, A&A, 507, 355Falgarone, E., & McKee, C.F. 2010, in preparationFatuzzo, M. & Adams, F.C. 2002, ApJ, 570. 210Fiedler, R. A. & Mouschovias, T. Ch. 1992, ApJ, 391, 199Fiedler, R. A. & Mouschovias, T. Ch. 1993, ApJ, 415, 680Gammie, C. F., Lin, Y. T., Stone, J. M., & Ostriker, E. C. 2003,ApJ, 592, 203Goodman, A. A., Barranco, J. A., Wilner, D. J., & Heyer, M. H.1998, ApJ, 504, 223Heiles, C. & Troland, T. H. 2005, ApJ, 624, 773Hennebelle, P., & Chabrier, G. 2008, ApJ, 684, 395Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25Heyer, M. H. & Brunt, C. M. 2004, ApJ, 615, 45Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ,699, 1092Iroshnikov, P. S. 1963, AZh, 40, 742 (English transl. SovietAstron., 7, 566 [1964])Klessen, R. S., Heitsch, F., & Mac Low, M.-M. 2000, ApJ, 535,887Klessen, R. S., & Hennebelle, P. 2009, arXiv:0912.0288Kraichnan, R. H. 1965, Phys. Fluids, 8, 1385Krumholz, M. R., McKee, C. F. 2005, 630, 250Kulsrud, R. & Pearce, W. P. 1969, ApJ, 156, 445Larson, R. B. 1981, MNRAS, 194, 809Li, P. S., Norman, M. L., Mac Low, M.-M., & Heitsch, F. 2004,ApJ, 605, 818Li, P. S., McKee, C. F., & Klein, R. I. 2006, ApJ, 653, 1280(LMK)Li, P. S., McKee, C. F., Klein, R. I., & Fisher, R. T. 2008, ApJ,684, 380 (LMKF)Lizano, S. & Shu, F. H. (1989), ApJ, 342, 834Mac Low, M.-M., Norman, M. L., Konigl A., & Wardle, M. 1995,ApJ, 442, 726 Mac Low, M.-M. & Smith, M. D. 1997, ApJ, 491, 596Mac Low, M.-M. 1999, 524, 169Mac Low, M.-M., & Klessen, R. S. 2004, Reviews of ModernPhysics, 76, 125Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364McKee, C. F. 1989, ApJ, 345, 782McKee, C. F., & Ostriker, E. C. 2007, ARAA, in press.Mestel, L., & Spitzer, L. 1956, MNRAS, 116, 503Mouschovias, T. Ch. 1976, ApJ, 207, 141Mouschovias, T. Ch. 1977, ApJ, 211, 147Mouschovias, T. Ch. 1979, ApJ, 228, 475Mouschovias, T. C. 1987, NATO ASIC Proc. 210: PhysicalProcesses in Interstellar Clouds, 453Mouschovias, T. C., & Psaltis, D. 1995, ApJ, 444, L105Mouschovias, T. C., & Spitzer, L., Jr. 1976, ApJ, 210, 326Myers, P. C., & Khersonsky, V. K. 1995, ApJ, 442, 186Myers, P. C., & Lazarian, A. 1998, ApJ, 507, L157Nakano, T., & Nakamura, T. 1978, PASJ, 30, 671Nakano, T. & Tademaru, E. 1972, ApJ, 173, 87Nakamura, F. & Li, Z. Y. 2008, ApJ, 687, 354Onishi, T., Mizuno, A., Kawamura, A., Tachihara, K., & Fukui,Y. 2002, ApJ, 575, 950Ostriker, E. C., Gammie, C. F., & Stone, J. M. 1999, ApJ, 512,259Padoan, P. 1995, MNRAS, 277, 377Padoan, P., & Nordlund, ˚A. 1999, ApJ, 526, 279Padoan, P. & Nordlund, ˚A. 2002, ApJ, 576, 870Padoan, P., Willacy, K., Langer, W., & Juvela, M., 2004, ApJ,614, 203Padoan, P., Nordlund, ˚A., Kritsuk, A. G., Norman, M. L., & Li,P. S. 2007, ApJ, 661, 972Pan, L. & Padoan, P. 2008, submitted, astro-ph 0806.4970Passot, T., V´azquez-Semadeni, E.,& Pouquet, A. 1995, ApJ, 455,536Plume, R., Jaffe, D. T., Evans, N. J., II, Martin-Pintado, J., &Gomez-Gonzalez, J. 1997, ApJ, 476, 730Shu, F. H. 1983, ApJ, 273, 202Solomon, P. M., Rivolo, A. R., Barret, J. & Yahil, A. 1987, ApJ,319, 730Spitzer, L., Jr. 1968, Diffuse Matter in Space (New York:Interscience)Tachihara, K., Onishi, T., Mizuno, A., & Fukui, Y. 2002, A&A,385, 909Tassis, K. & Mouschovias, T. Ch. 2005, ApJ, 618, 769Tassis, K. & Mouschovias, T. Ch. 2007, ApJ, 660, 388Tilley, D. A., & Pudritz, R. E. 2004, MNRAS, 353, 769Tilley, D. A. & Pudritz, R. E. 2007, MNRAS, 382, 73.Tomisaka, K., Ikeuchi, S., & Nakamura, T. 1988, ApJ, 335, 239T´ o th, G. 1995, MNRAS, 274, 1002Troland, T. H., & Crutcher, R. M. 2008, ApJ, 680, 457V´azquez-Semadeni, E., Ballesteros-Paredes, J., & Klessen, R. S.2003, ApJ, 585, L131V´azquez-Semadeni, E., Kim, J., Shadmehri, M., &Ballesteros-Paredes, J. 2005a, ApJ, 618, 344 V´azquez-Semadeni, E., Kim, J., & Ballesteros-Paredes, J. 2005b,ApJ, 630, L49V´azquez-Semadeni, E., Gonz´alez, R. F., Ballesteros-Paredes, J.,Gazol, A., & Kim, J. 2008, MNRAS, 390, 769Wakelam, V., & Herbst, E. 2008, ApJ, 680, 371Ward-Thompson, D., Andr´e, R., Crutcher, R., Johnstone, D.,Onishi, T., WIlson, C. 2007, in Protostars and Planets V, ed.Reipurth, B., Jewitt, D., & Keil, K., p.33Williams, J. P., De Geus, E. J., & Blitz, L. 1994, ApJ, 428, 693Williams, J. P., Bergin, E. A., Caselli, P., Myers, P. C., & Plume,R. 1998, ApJ, 503, 689 Williams, J. P., Blitz, L., & McKee, C. F. 2000, in Protostars andPlanets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell, p.97Zuckerman, B., & Evans, N. J. 1974, ApJ, 192, L149Zuckerman, B., & Palmer, P. 1974, ARA&A, 12, 279Zweibel, E. G. 2002, ApJ, 567, 962Zweibel, E. G. & Brandenburg, A. 1997, ApJ, 478, 563
HD Turbulence Simulations with Ambipolar Diffusion 25 − − − − − − − − M A , i R AD,c
Fig. 1.—
Clumps’ AD Reynolds number, R AD , c , versus ion Alfv´en Mach number squared, M A i , for models m3c2r-1(circles),m3c2r1(squares), and m3c2r3(diamonds) at the end of the simulation. The straight line shows R AD , c = M A i and the solid symbolsindicate the values of R AD and M A i for the whole box. − . − − . − − . − − . . . . l og M c / M log (d N / d ln M c /M ) . . . l og k c k i n , m a x k i n , m a x k s Γ f it = − . ± . Fig. 2.—
Clump mass functions for AD models of R AD = 1200 with grid sizes of 256 (dashed line) and 512 (m3c2r3, solid line). Thesonic wave number, k s , and the minimum scale of the inertial range, k in , max , for model m3c2r3 are plotted as vertical solid lines; k in , max for the 256 model is plotted as a vertical dashed line. The clump wavenumber k c based on equation (40) is plotted at the top of the figurefor reference; note that k increases to the left. See § HD Turbulence Simulations with Ambipolar Diffusion 27
Fig. 3.—
3D spatial distribution of clumps, identified by CLUMPFIND with minimum mean radius of 6 cells, from model m3c2r1.Different gray-scale shadings (different colors in the online version of the paper) are used only to visually separate overlapping individualclumps. − . − − . − − . − . . . . l og M c / M log (d N / d ln M c /M ) Γ f it = − . ± . Γ f it = − . ± . Fig. 4.—
Clump mass functions for models m3c2r-1 [ R AD ( ℓ ) = 0 .
12; solid line], m3c2r1 [ R AD ( ℓ ) = 12], and m3i (ideal MHD). Thedashed line and dot-dashed line show the best fitting higher-mass slope Γ fit for models m3c2r1 and m3i, respectively. Model m3c2r-1, whichhas the strongest ambipolar diffusion, has the steepest higher-mass slope (see § R AD ( ℓ ) = 1200], which is not shown. All clumps with radius larger than 3 cells are included in the plot. HD Turbulence Simulations with Ambipolar Diffusion 29 − − − M c / M ( µ Φ ,c / µ Φ ,0 ) / (M c /M ) Fig. 5.—
Convergence study for the mass-to-flux ratio of the clumps. Values of normalized mass-to-flux ratios of clumps in the R AD = 1200model are plotted against normalized clump mass, for resolutions of 256 (blue cirlces) and 512 (red squares; model m3c2r3). By plotting( µ Φ ,c /µ Φ , )/( M c /M ) / versus M c /M , data points are projected horizontally for easy visual comparison. The straight line shows themean of the normalized mass-to-flux ratio for elliptical clumps using the mean density and B-field of the clumps from model m3c2r3. See § . . − − − M c / M ( µ Φ ,c / µ Φ ,0 ) / (M c /M ) Fig. 6.—
Normalized mass-to-flux ratios, ( µ Φ ,c /µ Φ , )/( M c /M ) / , of clumps in models m3c2r-1 (blue circles), m3c2r1 (black crosses),and m3c2r3 (red triangles) plotted versus normalized clump mass M c /M . The mean values of h µ Φ ,c /µ Φ , / ( M c /M ) / i for the threemodels are plotted as the horizontal lines ( R AD ( ℓ ) = 0 .
12 blue solid, R AD ( ℓ ) = 12 black dashed, and R AD ( ℓ ) = 1200 red dot-dashed).The model with the largest value of R AD has a slightly lower average mass-to-flux ratio and a smaller dispersion of the mass-to-flux ratios. HD Turbulence Simulations with Ambipolar Diffusion 31 R AD,c r c /l −5 −5 −3 −1 R AD,c ρ i , c / ρ −5 R AD,c ρ n , c / ρ −5 R AD,c U B , c / U B , −5 −5 −4 −3 −2 R AD,c M c / M −5 −2 R AD,c χ i , c / χ i , (c) (b)(e) (f)(d)(a) Fig. 7.—
Other normalized physical properties of clumps as functions of the AD Reynolds number of the clumps, R AD , c , for the modelsm3c2r-1 (on the left) and m3c2r3 (on the right): (a) radius, r c /ℓ , (b) ion density; h ρ i i c / ¯ ρ ; (c) neutral density, h ρ n i c / ¯ ρ ; (d) magneticenergy density, U B,c /U B, ; (e) clump mass, M c /M ; and (f) ionization mass fraction, χ i,c /χ i, . See § TABLE 1AD Reynolds Number R AD for Observed Molecular Clumps (Crutcher 1999) Cloud β log n R M M A T k R a AD (H cm − ) (pc) (K)W3 OH 0.07 6.8 0.02 1.9 0.3 100 3.0DR 21 OH1 0.21 6.3 0.05 4 1.3 50 37.3Sgr B2 0.0008 3.4 22 22 0.4 70 10.3M17 SW 0.008 4.5 1 7 0.5 50 6.3W3 (main) 0.13 5.5 0.12 4.8 1.2 60 24.1S106 0.04 5.3 0.07 3.6 0.5 30 3.7DR 21 OH2 0.41 6 0.05 4 1.8 50 51.5OMC-1 0.65 5.9 0.05 1.7 1 100 21.9NGC 2024 0.35 5 0.2 3.7 1.6 25 72.7S88 B 0.056 3.8 0.7 5.9 1 40 12.9B1 0.17 4 0.2 3.6 1.1 12 15.7W49 B 0.024 3 1 5.9 0.6 10 6.3W22 0.033 3 4 3.5 0.5 10 20.5W40 0.027 2.7 5 10 1.2 10 42.4 ρ Oph 1 0.42 3.2 0.8 3.5 1.6 25 41.6OMCN-4 > > > > > > > > > > > > ρ Oph 2 > > > > > > > > > > > > > > > > > > > > > > > > a R AD computed using equation (15) HD Turbulence Simulations with Ambipolar Diffusion 33
TABLE 2Model Parameters and Regimes of AD
Model a γ AD R AD ( ℓ ) R AD ( ℓ ) bt Regime of ADm3c2r-1 4 0.12 0.076 IIIm3c2r0 40 1.2 0.70 II ∼ IIIm3c2r1 400 12 10.1 IIm3c2r2 4000 120 103.2 IIm3c2r3 40000 1200 1022 Im3i ∞ ∞ I a Models are labeled as “mxcyrn,” where x is the thermal Mach number, y = | log χ i | , and n = log( R AD ( ℓ ) / . b R AD from models using time-dependent ionization (see § c Root mean squared (rms) values.
TABLE 3Comparison of Clump Properties in Models with Different R AD Model m3c2r-1 m3c2r1 m3c2r3 R AD ( ℓ ) 0.12 12 1200 n vn ( k ) a . ± .
02 1 . ± .
03 1 . ± . b fit - − . ± . − . ± . h µ Φ ,c i /µ c Φ , . ± .
004 0 . ± .
005 0 . ± . σ ( h µ Φ ,c i /µ Φ , ) d .
056 0 .
074 0 . h ρ c i . ± .
24 5 . ± .
23 3 . ± . h R cz /R c, ⊥ i . ± .
02 1 . ± .
02 2 . ± . h r c i (cell) e . ± . ± . . ± . h N c i ( d c >
12 cells) f
698 434 349 a Velocity power spectral index of neutral component. b Slope of the ClMF in the inertial range. The data for model m3c2r-1 do not have a single power law over the inertial range. c Clump mass-to-flux ratio normalized by that of the whole box. d Dispersion of clump mass-to-flux ratio normalized by that of the whole box. e Mean radius of clumps in units of number of cells. ff