accepted toApJ January4th,2014 PreprinttypesetusingLATEXstyleemulateapjv.04/17/13 BONDI-HOYLE ACCRETION IN AN ISOTHERMAL MAGNETIZED PLASMA Aaron T. Lee1, Andrew J. Cunningham2, Christopher F. McKee1,3, Richard I. Klein1,2 accepted to ApJ January 4th, 2014 ABSTRACT Inregionsofstarformation,protostarsandnewbornstarswillaccretemassfromtheirnatalclouds. Thesecloudsarethreadedbymagneticfieldswithastrengthcharacterizedbytheplasmaβ—theratio 4 ofthermalandmagneticpressures. Observationsshowthatmolecularcloudshaveβ < 1,somagnetic 1 fields have the potential to play a significant role in the accretion process. We ha∼ve carried out a 0 numerical study of the effect of large-scale magnetic fields on the rate of accretion onto a uniformly 2 moving point particle from a uniform, non-self-gravitating, isothermal gas. We consider gas moving n with sonic Mach numbers of up 45, magnetic fields that are either parallel, perpendicular, or a oriented 45◦ to the flow, and β aMs lo≈w as 0.01. Our simulations utilize adaptive mesh refinement in J order to obtain high spatial resolution where it is needed; this also allows the boundaries to be far 7 from the accreting object to avoid unphysical effects arising from boundary conditions. Additionally, 2 we show our results are independent of our exact prescription for accreting mass in the sink particle. Wegivesimpleexpressionsforthesteady-stateaccretionrateasafunctionofβ and fortheparallel ] M A and perpendicular orientations. Using typical molecular cloud values of 5 and β 0.04 from the literature, our fits suggest a 0.4 M star accretes 4 10−9 M /yMear,∼almost a fa∼ctor of two ⊙ ⊙ G ∼ × less than accretion rates predicted by hydrodynamic models. This disparity can grow to orders of . magnitude for stronger fields and lower Mach numbers. We also discuss the applicability of these h accretionratesversusaccretionratesexpectedfromgravitationalcollapse,andunderwhatconditions p a steady state is possible. The reduction in the accretion rate in a magnetized medium leads to an - o increase in the time required to form stars in competitive accretion models, making such models less r efficient than predictedby Bondi-Hoylerates. Our results shouldfind applicationin numericalcodes, t s enabling accurate subgrid models of sink particles accreting from magnetized media. a Subject headings: ISM: magnetic fields — magnetohydrodynamics (MHD) — stars: formation [ 1 v 1. INTRODUCTION induces global gravitational collapse, resulting in super- 0 Accretionisubiquitousinastrophysics. Withexamples sonicinfalleitherdirectlyontothestellarsurfaceorintoa 1 including protostellar accretion from molecular clouds, surrounding centrifugally supported disk. Material that 0 masstransferbetweenbinarycompanions,andgasfalling ultimately ends up on the star comes from a local gravi- 7 onto a supermassive black hole in the center of galac- tationallyboundregionoftheparentmolecularcloud. If . the core is not collapsing directly onto the star+disk, or 1 tic nuclei, understanding how (or whether) a gravitating if the core is exhausted and the star is moving through 0 source gathers mass has received much attention over the more tenuous regions of the cloud, another accre- 4 the pastcentury. Inthe caseofstarformation,consider- tion mechanism is at play. Here the local gas initially 1 able study has been given to understanding the process unbound to the star can be captured and subsequently : of accretion from a background medium. Knowing how v accreted. The self-gravity of this local gas is negligible much mass a star can accrete from its natal cloud will i relative to the gravity of the star itself. Such accretion X help elucidate, for example, whether the final mass of is often called Bondi accretion when the star is station- the star is determined primarily through gravitational r ary or Bondi-Hoyle(-Lyttleton) accretion when the star a collapse (e.g., Shu 1977) or through post-collapse accre- ismovingrelativetothebackgroundgas,namedafterthe tion (e.g., Bonnell et al. 1997, 2001). Mass accretion pioneering investigators(Hoyle & Lyttleton 1939;Bondi also could play a role in the dynamics of stars in clus- & Hoyle 1944; Bondi 1952). ters. Iftheaccretorismovingrelativetothebackground Theprimarygoalofthisworkisconcernedwithunder- gas, then the accretion of mass and momentum will be standing the steady-state mass accretionrate for Bondi- non-spherical, and this may play a role in the radial re- Hoyle accretion when the background gas is an isother- distribution of objects in stellar clusters (Lee & Stahler mal plasma pervaded by a magnetic field. In particular, 2011). we seek to construct an interpolation formula that re- Several physical processes exist for transferring mass produces both known analytic and numerical results as from the cloud to the surface of a (proto)star. In core- wellasthesteady-stateaccretionrateswewillobtainvia collapse models (Shu 1977), a dense core’s self-gravity numerical simulations. In our work and these previous a.t.lee@berkeley.edu works, the effects of stellar winds and outflows are ne- 1DepartmentofAstronomy,UniversityofCaliforniaBerkeley, glected. We begin this study by summarizing some of Berkeley,CA94720 2Lawrence LivermoreNational Laboratory, P.O.Box808, L- the known results in the next section. From there, we 23,Livermore,CA94550 propose new interpolation formulas for the mass accre- 3Department of Physics, University of California Berkeley, tion rate of magnetized Bondi-Hoyle flow. This function Berkeley,CA94720 2 Lee et al will have two free parameters, which we fix by fitting to where we have introduced the sonic Mach number M≡ numericalsimulations. Section3discussesourmethodol- v /c . The characteristicvelocityfor Bondi-Hoyleaccre- 0 s ogy and numerical convergence studies of the numerical tion is code. Section 4 presents the numerical results and the v =(c2+v2)1/2 , (7) BH s 0 results of fitting our proposed interpolation formulas to and the corresponding Bondi-Hoyle radius is the simulationdata. InSection5wediscussthe applicability of such steady-state models in regions of star for- GM r ∗ B mation. Section 6 concludes this work with a summary rBH = v2 = 1+ 2 . (8) and discussion. How these models can be implemented BH M in sub-grid and sink particle algorithms is discussed in We will see in Section 4 that magnetic fields reduce the Appendix A. accretion rate below these values. Furthermore, for the fiducial values of n and T and for M > 0.4M , Bondi 0 ∗ ⊙ 2. MASSACCRETIONRATES accretion is not in a steady state (Section 5). 2.1. Known Results We shall express all accretion rates in terms of the Bondiaccretionrate in two equivalentforms. For exam- The study of steady-state accretion from an initially ple, the Bondi-Hoyle accretion rate will be written as uniform background medium has enjoyed many analyt- ical and numerical studies. Edgar (2004) gives a nice 3 c pedagogical review of some of the earlier work. Hoyle M˙ =φ 4πλr2 ρ v =φ s M˙ . (9) BH BH· BH 0 BH BH v B & Lyttleton (1939) first solved the problem for a point (cid:18) BH(cid:19) particle of mass M moving through a collisionless fluid ∗ This first form emphasizes the underlying physical pa- at(hypersonic)speedv . Matterwasfocusedintoavan- 0 rameters, and we have introduced a correction factor ishinglythinwakeandaccretedthroughaspindle down- φ =φ ( ),whichwillbeoforderunity. Thesecond BH BH stream of the accretor. The accretion rate was M form is 3 M˙HL =4πrH2Lρ0v0 = 4πG2vM03 ∗2ρ0 , (1) M˙BH =(cid:18)vBcHs,eff(cid:19) M˙B . (10) for the far-field mass density ρ0. Associated with v0 is Here, the effective Bondi-Hoyle velocity vBH,eff is an in- the characteristic radius terpolationformula;therationaleforintroducingthesec- ond form will become clear below. GM r ∗ , (2) SimulationshaveshownthatBondi’sinterpolationfor- HL ≡ v02 mula(φBH =1/λ)canbeinerrorbyseveralten’sofper- cent (Shima et al. 1985; Ruffert 1994). These authors, which measures the dynamic length scale within which among others, have considered the non-isothermal case gravity wins over the inertia of the gas. In the opposite aswellandhaveproposedtwo-dimensionalinterpolation limit ofstationaryor subsonicmotion, the thermalpres- formulas(in andγ)tomatchsimulationresults. Typ- sure exceeds the ram pressure of the gas by a factor of M ically such formulas are monotonically decreasing func- (c /v )2 for sound speed c . Bondi (1952)analytically ∼ s 0 s tions of both and γ and agree well the simulations. solvedthe problemforastationaryaccretor,arrivingat M A complication is that Ruffert (1994, 1996) has shown 4πλG2M2ρ thataccretionratesdonotdecreasemonotonicallyas M˙B =4πλrB2 ρ0cs = c3s ∗ 0 , (3) ianscyrmeapsteost,ebtuotMi˙nste.adFoinrcrtehaeseisontehaerrmMal∼cas1e,awnde hthaMevne HL which becomes found that M 2 n T −3/2 M (1+ 2)3/2[1+( /λ)2]1/2 =1.02×10−6 0.4 M∗ 104 c0m−3 10 K yr⊙(4.) φBH = M 1+ 4M , (11) (cid:18) ⊙(cid:19) (cid:16) (cid:17)(cid:18) (cid:19) M Here we have defined the Bondi radius corresponding to rB ≡ GcM2s ∗ =9.0×1016(cid:18)0.4MM∗⊙(cid:19)(cid:18)10TK(cid:19)−1 cm . MBH ≡ vBcHs,eff = [1(+1+(MM/λ4))12/]31/6 , (12) (5) and agrees with the numerical results of Ruffert (1996, In our numerical evaluations, we have normalized the forγ =1.01)andthosereportedbelow withamaximum temperature T to 10 K, the number density n to 104 0 error of 27%. Observe that φ 1 as 0 and particles per cm3, and masses to the solar mass M⊙. φ 1/λ for 1. ThBisHf→unction Mis p→lotted in The mass density is related to the number density by BH → M ≫ Figure 1. ρ =(2.34 10−24 grams) n . Thesymbolλisafunction 0 × · 0 Other numerical studies of accretion have studied of the adiabatic index γ (λ = exp(3/2)/4 1.12 for an ≈ the role of additional physics like radiation pressure isothermal gas, γ =1). (Milosavljević et al. 2009), turbulence (Krumholz et al. Both limits then established, Bondi proposed his ven- 2006),turbulenceandmagneticfieldsinsphericallysym- erable Bondi-Hoyle interpolation formula that connects metric accretion (Shcherbakov 2008), the presence of a the stationary and hypersonic regimes: disk (Moeckel & Throop 2009), or thermal instabilities 4πρ r2c M˙ /λ (Gaspari et al. 2013). Our finding is that Equation (9) M˙ = 0 B s = B , (6) withφ givenbyEquation(11)isareasonablemeasure BH (1+ 2)3/2 (1+ 2)3/2 BH M M Magnetic Bondi-Hoyle Accretion 3 Cunninghametal.(2012)expressedtheaccretionrateas β=0.01 β 1/2 β=0.1 M˙ = M˙ , (17) φBH φABH,⊥ β=1 A (cid:18)βch(cid:19) B 100 where β is a numerical factor; this corresponds to ch φ = (2/β )1/2. They estimated β 5, so that φ φABH,∥ φA =0.63. Fchromhere, they generalizedchthi≈s to include a A 100 φABH(M=0,β) finite temperature (the “Alfvén-Bondi” case). By writ- ing v , r v , r in Equation (15), the accretion A A AB AB → rate becomes 0.2 −6 −3 0 3 6 c 10-1 log10(β) M˙AB=φAB 4πλrBrABρ0vAB =φAB s M˙B(,18) 10-1 100 101 102 103 · (cid:18)vAB(cid:19) Mach Number c M = s M˙ , (19) Fig.1.—VariousnumericalparametersφasafunctionofMach v B number andplasma β. Colorsrefer to different magnitudes of (cid:18) AB,eff(cid:19) 1β/,λw.hTerheeMassutbhpelloitnesshtoywlesstdhieffestraftoirondaiffryerleimntitφo.fTφhAeBhHigahsMafulinmctiitoins where vAB ≡ (c2s +vA2)1/2 and rAB = GM∗/vA2B.4 The of β. This has the same functional form as φAB, defined in the effective Alfvén-Bondivelocity, vAB,eff, can be chosen to text, but with different fitting parameters βch and n, which are provideaninterpolationformulabetweentheAlfvénand determinedin 4.2. Theasymptotesarep2/βch 0.32andunity. Bondi cases that agrees best with the numerical simula- § ≈ tions; Cunningham et al. (2012) adopted n/2 1/n n/2 1/n of the accretion rate for isothermal Bondi-Hoyle accre- v cn+ βch vn = 1+ βch c , tion when additional physics do not play an appreciable AB,eff ≡" s (cid:18) 2 (cid:19) A# " (cid:18) β (cid:19) # s role in the dynamics of the gas near the accretor. (20) One physical effect that could play an important role which gives in the gas dynamics is a global magnetic field. In star- forming regions, there is ample evidence that molecular n/2 −1/n β clouds are pervaded by magnetic fields (Crutcher 1999; M˙ = 1+ ch M˙ . (21) AB B McKee & Ostriker 2007), whose strength (i.e., its abil- " (cid:18) β (cid:19) # ity to influence dynamics) can be characterized by the They found that n=0.42 and β =5.0 gave agreement plasma β, the ratio of the thermal pressure to the mag- ch with their numerical results to within 5% for β 0.01. netic pressure: ≥ ρc2 c 2 2 2.2. Alfvén-Bondi-Hoyle Accretion β s =2 s =2 MA , (13) We wish to extend the work of Cunningham et al. ≡ B2/8π v (cid:18) A(cid:19) (cid:18) M (cid:19) (2012)to the caseinwhichthe accretingmassis moving formagneticfieldamplitudeB. Wehavealsointroduced through a magnetized ambient medium. Our primary the Alfvén Mach number = v /v , the ratio of the interest is in star-forming regions, which are approxi- A 0 A gas velocity to the AlfvénMvelocity, v = B/√4πρ. Ob- mately isothermal because the dust and the molecules A servationsofthe Zeemaneffect, linearpolarizationemis- can efficiently radiate the energy supplied by compression of dust, and the Chandrasekhar-Fermimethod (see sion; we therefore assume that the gas is isothermal.5 the review of Crutcher 2012) have suggested molecular The magnetic flux in stars is orders of magnitude less clouds have β values of most order unity, but more have than that in the gas from which they form, so most of β <0.1 (e.g., Crutcher 1999). of the magnetic flux in the accreting gas decouples from Cunningham et al. (2012)have studied accretion from the gas and accumulates in the vicinity of the protostar a magnetized, isothermal, static medium. For the case (Zhao et al. 2011). As a result, even in cases where the in which thermal pressure is negligible (the low-β limit) thermal pressure ( ρc2) or ram pressure ( ρv2) ini- ∼ s ∼ they argued that gas would collapse along the field lines tially control the dynamics of the gas near the accretor, from a distance r above and below the point mass, accretioncanredistributemagneticfluxsothatthemag- B and would then a∼ccrete from an Alfvén radius, neticpressure( B2)eventuallydominatesthedynamics ∼ near the accret∝ing object. In steady-state Bondi accre- r GM∗ = c2s r (14) tion from a magnetized gas, Cunningham et al. (2012) A ≡ v2 v2 B found that even if β was initially > 1, a steady-state A A was reached when the gas within r of the accre- AB at velocity v . In our notation, tor had β 1. In a steady-state∼flow where there is A ∼ ≈ relative motion between the gas and the accretor (i.e., M˙A =φA 4πλrBrAρ0vA . (15) magnetized Bondi-Hoyle accretion), we anticipate that · As a result, the accretion rate varies as v−1 β1/2, 4 Note that Equation (16) of Cunningham et al. (2012) has a A ∝ typo,onefactorofrAB shouldbewrittenasrB instead. c β 1/2 5AswediscussinSection6,ourresultsshouldalsobeapplicable M˙ =φ s M˙ =φ M˙ . (16) tothecentral regionsofactivegalacticnuclei. A A B A B v 2 (cid:18) A(cid:19) (cid:18) (cid:19) 4 Lee et al if min(1, ), the inertia of the gas will play a inEquation(25). Onecanreadilyverifythatthishasthe A M ≪ M small role in setting the steady-state accretion rate, so correctlimits for Bondi, Bondi-Hoyle,and Alfvén-Bondi M˙ will be wellapproximatedby the Alfvén-Bondiresult accretion. We then obtain (Eq. 18). If instead 1, the inertia of the gas is MA ≫ n/2 −1/n able to drag away most of the magnetic flux so that the 1 1 β M˙ =min , n + ch M˙ . tchreetimonagrnaetteicshfioeuldldisapnportodaocmhitnhaenntoann-ymwahgenreetiazneddtBhoenadci-- ⊥ M3BH MBH "MBH (cid:18) β (cid:19) #  B Hoylelimit. Wewishtodevelopanapproximateanalytic (30) expressionfortherateofaccretionbyapointmassmov- An immediate interesting result of this formulation— ing at a constant speed through a uniform, isothermal, that is born out in our simulations, see 4.2—is for the § magnetized medium by further generalizing the above particular case of highly supersonic flow with an Alfvén knownresults. Ourexpressionwillthereforealsoinclude Machnumber > 1, the accretionratefor the perpendic- parameters n and β , which we can then adjust to best ular case reduce∼s to ch reproduce the results of our simulations of magnetized M˙ /λ Bondi-Hoyle flow. M˙ = B =M˙ ( 1, 1) , (31) WegeneralizeEquation(18)byreplacingr withr , ⊥ 3 HL M≫ MA ≥ B BH M rAB with rABH, and vAB with even when A 1. M ≈ v (c2+v2+v2)1/2 . (22) If the point mass is moving through a medium at an ABH ≡ s 0 A angle θ with respect to the magnetic field, we approxi- Here, and throughout the remainder of the paper, c is mate the accretion rate by s the isothermal sound speed. With M˙ M˙ cos2θ+M˙ sin2θ . (32) GM∗ ≃ k ⊥ r , (23) ABH ≡ v2 Indeed, we confirm for one of our simulations that when ABH θ = 45◦, the resulting accretion rate is decently approx- the accretion rate is then imated by the average of the two limiting rates. If the M˙ =φ 4πλr r ρ v orientation changes randomly in time, the proposed av- ABH ABH BH ABH 0 ABH · erage accretion rate is c3 =φ s M˙ (24) ABH(cid:18)vB2HvABH(cid:19) B M˙ ≃ 12 M˙k+M˙⊥ . (33) = v2 cv3s M˙B. (25) In order to test our pr(cid:16)oposed inte(cid:17)rpolation formulas, BH,eff ABH,eff! westudytheproblemofBondi-Hoyleaccretioninamag- netized plasma using the RAMSES MHD code (Teyssier Equations (24–25) do not take into account the orien- 2002)overarangeoffieldstrengthsandsonicMachnum- tation of the flow relative to the ambient magnetic field. bers relevant for star formation. These simulations em- InSection4.2below,wefindthatweneeddifferentinter- ploy the adaptive mesh refinement (AMR) capabilities polation formulas for the cases where the flow is parallel of the code to retain high spatial resolution where it is andperpendicularto themagnetic field. Forthe parallel needed—closetotheaccretingobject—whileallowingfor case, we generalize v to v with AB,eff ABH,k,eff alargecomputationaldomainto preventthe boundaries v v n β n/2 1/n ofthe domainfrominfluencing the steady-stateflow. As ABH,k,eff BH,eff ch noted above, we do not consider the effects of stellar = + . (26) cs "(cid:18) cs (cid:19) (cid:18) β (cid:19) # winds or outflows on the accretion rate. The numerical methodology is described in the next section. The accretion rate in this case is 1 β n/2 −1/n 3. NUMERICALMETHODS M˙ = n + ch M˙ , (27) Our methods are similar to those described in Cun- k M2BH (MBH (cid:18) β (cid:19) ) B ningham et al. (2012). Here we summarize the methods, highlighting the significant differences in this work, where was defined in Equation(12). Observe that MBH presenttheresultsofourconvergencestudy,andpointto this expression reduces to Bondi accretion for = β−1 0, to Bondi-Hoyle accretion if β , aMnd to wherethe readercanfindadditionaldetailsifinterested. → → ∞ WesolvetheequationsofidealMHDforanisothermal Alfvén-Bondi accretion for =0 and arbitrary β. The M gas with a fixed point mass at the origin, particularly factor φ in this case is then ABH,k mass conservation, v 2/3 ABH ∂ρ φABH,k =φBH · v ; (28) + ρv= SM(x) , (34) (cid:18) ABH,k,eff(cid:19) ∂t ∇· − it is also plotted in Figure 1. momentum conservation, For flows perpendicular to the field, we obtain better agreement with our simulations with the less-intuitive ∂ρv B2 + (ρvv)= P + interpolation ∂t ∇· −∇ th 8π (cid:18) (cid:19) v v v B B GM ρ ABH,⊥,eff max BH,eff, ABH,k,eff (29) + ·∇ ∗ xˆ S (x) v , (35) c ≡ c v 4π − x2 − M · s (cid:20) s BH,eff (cid:21) Magnetic Bondi-Hoyle Accretion 5 the induction equation, 2 ∂B (v B)=0 , (36) ∂t −∇× × 1 and the equation of state, ) P =ρc2 . (37) M0 th s ( 0 tHoertehevsisintkhepvaerltoiccliet,yaonfdthBegiass,thxeismtahgenpeotsicitifioenldr.elaSteivlfe- Log1−1 MA> 1A< 1 M gravity of the gas is neglected. In the code, the point −2 mass is represented by a fixed sink particle located at the center of the computational domain. The term S M −3 allows for mass accretion onto the central point mass if −4 −3 −2 −1 0 1 2 3 gasflows into a sphere of radius 4∆x, where ∆x is equal Log (1/β) 10 to the size of the grid cell on the finest AMR level. The accretedgas’s momentum is also removedfromthe grid, Fig.2.— Parameter space to be studied. Black dots represent models explored inthis work, with the two runs with arrows cor- thoughtheparticleisheldstationaryatthecenterofthe responding to β = 1030. Red dots are the stationary models of domain.6 The sink particle is allowed to accrete mass Cunninghametal.(2012). Withthischoiceofaxes,theleftverti- but not magnetic flux, and it accretes as much mass in cal axis approximates non-magnetic flow, where the bottom hori- zontal axis approximates stationary flow. The diagonal line plots a timestep ∆t as it can without introducing a new local A=1,whilethehorizontal lineplots =1. Ourrunsexplore maximum in the Alfvén speed amongst the cells located tMworegionsofthisparameterspacequitewMell. Theregionof <1 within a shell with radius r between 4∆x and 6∆x from and A > 1, not explored by us, was studied by RuffertM(1994, the accreting particle. That is, 1996)Min his investigations of non-magnetized isothermal Bondi- Hoyleaccretion. Since typical star formingregions have A 1 (Crutcher 2012), we also explore two cases with an Alfv´Men M≈ach 1 B max 0 ,ρ if x <4∆x numberofunity. SM(x)=∆t − 4πvA2,max! | | , thestellarfieldstrengthafewµGaussat 10AU,which 0 otherwise is already smaller than the field in the I∼SM. Therefore, (38) we neglect the field generated by the star itself. where v  = max(v (x); 4 x/∆x 6). The For all our integrations, the gas is initially uniform A,max A reader can also see the paragrap≤h c|on|taining≤Equation with density ρ0 and sound speed cs. The magnetic field (7) of Cunningham et al. (2012) for more details on the is initially set to be uniform in the zˆ direction with a − sink particle algorithm. magnitude set by β. The speed of the gas is initially set Since the sink particle accretes mass but not flux, the to v0, which is oriented either parallel or perpendicular cellsinteriortothesinkparticleradiusdecouplesthegas to the B-field, except in one case where we orient v at form the field. In reality, non-ideal MHD effects remove an angle of 45 degrees. We explore a parameter space of themajorityoftheaccretedgasfromthefieldwithinthe β and . We considerβ valuesof1030, 102, 10, 1, 0.1, M accretiondisk<100AUfromthestar;seeequation(48) and0.01andsonic Machnumbers thatrangefrom0.014 of (Li & McKee 1996)or the review of Armitage (2011). to 44.7. For a given β, we select our velocities to be Our sink particles will typically have a radius of 500 either equal to the Alfvén velocity or -1, -1/2, 1/2, or AU,soourtreatmentofnon-idealMHDeffectsrequ∼iresa 1 decade from this value. This gives us a combination sub-grid model; our prescription was given above.7 Fur- of runs that are both sub and super-sonic as well as thermore, this also means that gas just interior to the sub and super-Alfvénic. Table 1 tabulates the param- sink particle radius could be artificially affected by non- eter space explored, and Figure 2 shows this parameter idealeffects. Nonetheless,boththeexactprescriptionfor space graphically. The plot identifies four regions of pa- howgasisremovedfromthefieldlinesandthesizeofthe rameter space, depending on whether and A are M M sink particleareunimportantaslongas the gasentering greaterorlessthanone. Includingthestationaryrunsof thesinkregionhasacceleratedtofree-fall. Ifthishasoc- Cunningham et al. (2012), our runs explore two of these curred, the in-falling gas has causally disconnected from regions quite well ( , A > 1 and , A < 1). The M M M M the surrounding medium and any artificial prescriptions emptyregion( <1, A >1)isexploredbyhydrody- M M cannotalterthe far-fieldgas. We discusshowareresults namic models of Bondi-Hoyle flow (Ruffert 1996). The are independent of the sink particle conditions in more final region ( > 1, A < 1) is only explored by one M M detail in Section 4.2. simulation. Typical star forming regions have A 1 M ≈ In addition to the very small magnetic flux the star (Crutcher 2012), so we also explore two cases with an gains by accretion, the star might also generate its own Alfvén Mach number of unity. fieldthroughdynamoaction. Thefieldsofnewbornstars We carried out our computations using the RAMSES are observed with strengths of order kGauss, but the code (Teyssier 2002), an adaptive-mesh-refinement dipolecomponentofthefieldfallsoffas(R /R)3,making (AMR) code with an oct-tree data structure. The com- ∗ putational domain is a three-dimensional Cartesian do- tio6nsWoefnCoutnentihnagthathmeeatbsaeln.c(e20o1f2t)hies−aStyMp(oxg)r·avphtiecraml ministthaekee.qua- mcraetinizewitthheadolemnagitnh oonft5o0arBCianrteeasciahndbiraescet-iloenv.elWgreiddios-f 7Non-idealeffectscanalsoplayaroleatlargerdistance(.1000 643. Denote this levelaslevel L=0. We allow for seven AU)withinshocksthatoriginatefromthecollisionofin-fallinggas andthemagneticfieldthathasbeenfreedfromaccreted material additionallevelsofrefinement(L=1,2,...,7),with each (Li&McKee1996). level incrementing the grid-cell density by 23 above the 6 Lee et al TABLE 1 SimulationParameters β M MA rABH/∆x tend/tB hMfastia|| hMfastia⊥ (M˙/M˙B)|| (M˙/M˙B)⊥ (M˙/M˙B)45◦ 100 0.014 0.1 161 8 1.4 1.7 0.323 0.379 100 1.41 10 54 9 2.2 1.6 0.363 0.332 10 1.41 3.2 51 7 1.4 1.1 0.273 0.294 10 4.47 10 8 4 3.8 3.5 0.012 0.011 1 1.41 1 33 5 1.6 1.0 0.106 0.182 0.116 1 4.47 3.2 7 3 2.1 2.0 0.013 0.012 0.1 0.447 0.1 8 3 1.7 1.5 0.064 0.061 0.1b 4.47 1 16 0.5 1.8 1.1 0.00163 0.0112 0.1c 44.7 10 21 3 10−4 8.2 8.6 10−5 8.12 10−6 0.01d 4.47 0.32 2 ×0.5 1.3 1.00 0.0024 0.0×026 e 1.41 n/a 55 5 6.7 0.4 ∞e 4.47 n/a 8 5 5.5 0.01 ∞ a Computed as thevolumeaverage overthecells 5and 6 ∆x from thesink particle. Cells are included if the gas is in-falling (i.e., if v·x<0.) b Forthissimulation,twoadditionallevelsofrefinementareallowedfortheparallelrun,reducingthevalueof ∆xbya factor of 22. For theperpendicular run,only one additional levelis allowed, reducing∆x byafactor of 2. c For this simulation, one additional level of refinement is allowed and each dimension of the computational domain is reduced from 50rB to 50/64rB by a factor of 27, reducing thevalueof ∆x by a factor of 27+1. d Forthissimulation, oneadditional levelof refinement is allowed, reducingthevalueof ∆x byafactor of 2. e For thesesimulations, β is set to 1030 to approximate non-magnetic flow. previous level. Each grid cell in the domain is initially B refined up to level L if its distance x from the center of M ˙ the domain satisfies / β=1.0, =1.41 2L2+51 rB 1. Figure 7 plots the shock front location as a function of time for β =1 and =1.41. A M least-squares fit to an exponential function suggests the shockwillvanishat 3r ,where dropsto unity. ABH A ∼ M Asnotedabove,switch-onshocksoccuronlyforβ <2. Our simulations show that perpendicular field compo- nents can develop in the flow at a finite distance downstream of the shock; of course, if the shock is not ex- actly parallel, then the upstream perpendicular component of the field can be amplified by the compression in the shock. Ineithercase,the perpendicular fieldcompo- nent causes material to pile up on one side of the sink particle,andtheinertiaofthegasdriveskinksinthefield farther downstream from the shock (see the left side of Figure 4). As the field piles up, magnetic pressure even- tually dominates and the field straightens itself out, but overshoots, collecting on the opposite side of the accretor. Theresultingmorphologyupstreamisanoscillating motion of the field and the flow with a period of t . B ∼ Fig. 5.— Same as Figure 4 but for perpendicular orientations. Material not immediately upstream of the accretor ImageresolutionreducedforArXiv. flows into the Mach cone and then primarily travels throughtheregionimmediatelydownstreamoftheshock = 1.41 and β = 1 in Figure 5; here the downstream towards the sink particle. However, we note that for M fieldlines wererecentlyreleased). While the mass accre- these particular simulations, the region of reduced B tion rate reaches an approximate steady state, the mor- along the Mach cone is most likely due to numerical re- phologyoftheflowshowslargerfluctuations. Thedetails connection because the field flips orientation over a few of this morphology depend on the initial orientation of gridcells;wehavenotcarriedouthigherresolutionsimu- the field to the flow, which we now consider in turn. lationstotestifthisstructureisconverged. However,we haveensuredthatthemassaccretionrateisconverged,as 4.1.1. Parallel Orientations discussed in the next subsection. The Machcone results In the parallel case, there are two types of shock: hy- inthedownstreambeinglessdensethanthebackground, drodynamic, in which B is unaffected, and switch-on rather than material forming a downstream wake as in shocks, in which a perpendicular component of the field the hydrodynamic limit. appearsbehind the shock front(Draine & McKee 1993). For =4.47, the bowshock only forms for β =0.01. M The conditions for the occurrence of a switch-on shock Even though β = 0.1 gives = 1 initially, the re- A M are (1) > 1: (2) v > c , corresponding to β < 2; gion interior to r is so close to the accretor that any A A s ABH M and (3) the post-shock flow must be less than v cosθ , decrease with β also occurs with an increase in , re- A 2 M where the subscript 2 denotes post-shock quantities and sulting in < 1 never being satisfied. For larger β, A M θ is the angle between the magnetic field and the flow even less field enhancement occurs upstream. As a re- velocity. The first two requirements ensure that the sult, the = 4.47 runs resemble non-magnetic flows, M shock velocity exceeds the fast-wave velocity v , which with a downstream “wake” forming as the region of gas F is max(v , c ) = v in this case. The third, post-shock shocked from the Mach cone. The majority of the mass A s A requirement translates to an upper bound on the pre- accretion occurs through this wake. Magnetic Bondi-Hoyle Accretion 9 ) 1.0 5 B r β=1, =1.41 dius (0.8 M aP(O=enr−lyp0 .p.0 o1fii,ntbt==s aw0.+i2t2hb tt≥tB fitted) 4 a R0.6 3 k )H c B o Par. fit=a[1 bexp( t/c)] A Sh0.4 a=0.59,b=1.12,−c=1.55− 2(r m a e0.2 1 r t s p U 0.0 0 2 4 6 8 10 Time t/t B Fig.7.—Locationoftheshockfrontalongtheupstreamaxis. To providedataforafit,thelocationisplottedatt/tB=0.2,0.5,and thenevery0.5tBuntil3tB,whenthesimulationsend. Theparallel orientationisplottedwithcircles,perpendiculartriangles. Aleast- squaresfittothefunctionshownisperformed. Theparallelshock velocity tends to zero at 3rABH, where the flow is unchanged from the background flow∼. For the perpendicular case, the shock maintains a nearly constant speed equal to the local magnetofast velocity. flow, a weak shock moves outward at v . F ∼ The = 1.41 runs all show the development of a M dense irrotational disk around the accretor interior to 0.25r , with the disk normal perpendicular to the in- B ∼ coming flow. For β = 1 and 10, this disk attaches to a downstream wake. In the β = 100 case, colliding flows have made the inner 0.25r flow unstable, similar to B ∼ theoscillatingflowwesawintheparallelcases. Theweak fieldendsupdrapedaroundtheshocksthatformaround thesinkparticle. For =4.47,againtheflowresembles M the non-magnetic case. No bowshock is launched except in the β = 0.1 case, but, as discussed above, this is a transient of the flow. We also perform one run at β = 1, = 1.41 with M a 45-degree angle between the flow and magnetic field. Within r , the flow tends to align itself with the ABH ∼ local magnetic field, and the general flow resembles that of the parallel orientation. The remaining runs, with Fig.6.— Characteristic flow quantities for = 1.41,β = 1 at t = 3tB. The left and right columns shownMthe parallel and M1030=, 4a4re.7daynndaβmi=ca0ll.y1 odromwiinthatMed =by1.t4h1e, r4a.4m7 apnredssβur=e perpendicularorientations. Thetoprowshows ,thesecondrow shows β, the third row shows A, and the bMottom row shows of the gas and therefore closely resemble non-magnetic (B/B0)2. The colormaps are inMlog10 space, where the axes are hydrodynamic flow. linear(withunitsofrB). Theblackcircleindicatesthesizeofthe sinkparticle,equalto4∆x. ImageresolutionreducedforArXiv. 4.2. Mass Accretion Rate For each simulation, the mass accretion rate, M˙ , rises 4.1.2. Perpendicular Orientations withtimeandthenlevelsoffafterafewt =r /c . Each B B s For the perpendicular orientation, shocks occur if the simulationis rununtil the rate plateausfor a least a few flow velocity exceeds vF, which is (vA2 +c2s)1/2 = cs(1+ tABH. The ratequotedinTable 1is the averageoverthe 2/β)1/2 in this case. Even if this condition is not ini- last 1/6 of the integration, following Cunningham et al. tially satisfied, a magnetosonic wave launched from the (2012). sink particle boundary can steepen into a weak shock as Figure 8 gives an example of the time evolution of M˙ it moves upstream from the sink particle. We see this, for the β =1 runs. The rate is averaged over 0.2t bins B for example, in the = 1.41, β = 1 case. Initially to reduce noise. In general, when < 1 or 1, v = (3/2)1/2v . ImMmediately upstream of the shock, the orientation of the field makes liMttle differenMceAin≫the F 0 both ρ and B are increasedfrom compression,but ρ has final accretion rate despite the very different morpholo- increased even further from material falling down field gies. When 1, the perpendicular rate always ex- A lines. This results in an increased β and reduced vF. ceeds the paMrallel≈rate. When 1, M˙HL well ap- Figure 7 plots this front location. Since β has increased proximates the perpendicular acMcret≫ion rate, even, sur- to 1.2, the wave (which has steepened into a weak prisingly, when 1. sho∼ck)movesataspeed (1+2/1.2)1/2c v 0.22c In the high βM( A1≈00) subsonic runs of both our work s 0 s ∼ − ∼ ≥ upstream, relative to the accretor. Perpendicular to the and that of Cunningham et al. (2012), the accretor un- 10 Lee et al las generalized and built upon previous known analytic B M and numerical results. In particular, we have proposed ˙/ β=1.0 a simple interpolationformula for the Bondi-Hoyle limit M ˙ te 10-1 (REuqffuearttio(1n9996),).wWhiechpemrfaotrcmheeds ttwheo βsim=u1la0t3i0onrunressutoltsveor-f Ra =1.41 ify this equation, finding it underestimates the true ac- n M cretionratebyonly19%and3%for =1.41, 4.47,re- tio spectively. The massaccretionratefoMrthe diagonalcase ccre10-2 M=4.47 laiveesrbaegtewoefetnhtehpertehdeicptaedrapllaerlaalnledlpanerdppenerdpiceunldaircuralatresr;atthees A s Parallel Diagonal (Equation 33) reproduces this diagonal case to 18%. s Perpendicular We perform a least-squares fit to Equations (27) and a M10-3 (30) using the union of data from this work and from 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Cunninghamet al.(2012). Since the values ofthe accre- Time t/t B tionratecanvaryoverordersofmagnitude,wedefinethe Fig.8.— Mass accretion rates as a function of time and field residualsintheleast-squaresfunctiontobethedifference orientation for β = 1.0. All rates are normalized to the Bondi of the logarithms rather than of the absolute values: accretion rate (Equation 3). The short (red) lines identify the steady-state Bondi-Hoyleaccretionrates(Equation9). S = (log M˙ log M˙ )2 . (39) 10 data− 10 fit dergoes a period of rapid accretion before the accretion Each data poinXt is given an equal statistical weighting. rate suddenly drops to 1/2 of the original value (see MinimizingS,wefindβ =18.3 0.004andn=0.94 ∼ ch Figure 6, right panel of Cunningham et al. 2012). The ± ± 0.15withS =0.956. Wedonotincludethediagonalrun reason for this effect is that the magnetic field eventu- or the two hydro runs in our fit. allybecomesdynamicallydominantafterenoughfluxhas The standard errors show that matching the data to built up near the accretor. Cunningham et al. (2012) these interpolation formulas is not terribly sensitive to showed that this occurs after t > (β/100)1/2tB for a dy- the exact value of n. Since the data are consistent with namically weak field and could∼take an arbitrarily long using n = 1, we adopt this value for simplicity. Fixing time as β . We note that we do not see this effect n=1, a least-squares fit to the data yields β =19.8 for our →=∞1.41 run with β = 100. In this case, the 0.006 with S =0.96. We fix β as 19.8. ch ± M ch magnetic flux is unable to appreciably build up within In Section 2.2 we wrote all accretion rates in terms of rABH ∼ rB/4 before the gas pulls the flux downstream. M˙B,whichisconstantacrossourentireparameterspace. Below,whenwedeterminethebestfitparametersforour Since we are studying gas initially in uniform motion, interpolationformulas,weusetheinitial(larger)steady- normalizingtoM˙ (Equation10)is alsouseful. Indeed BH statevalueforthemassaccretionratesinceinastrophys- with this normalization, our parallel accretion rate can icalapplicationstheflowisoftennotsteadyforlongtime be written in terms of one parameter periods. Tworequirementsareneededtoensurethattheaccre- β1/2 =21/2 cs , (40) tion rate has converged. First, the value of M˙ should MBH MBH v (cid:18) A(cid:19) not depend on the resolution of the grid. As explained in 3, we have verified that this is the case. Increasing which varies as B−1 if the other parameters are held the§resolutionalsodecreasesthe sizeofoursink particle, constant ( BH is defined in Equation 12; recall that M whichhasaradiusof4∆x. Thesecondrequirementisto BH 1 as 0). Equations (27) and (30) can be M → M → ensurethatthesinkparticleboundaryconditionscannot rewritten as influence the valueofM˙ . Todo this,werequirethatthe M˙ 4.4 −1 k accreting gas pass through the fast magnetosonic point = 1+ , (41) at r > 4∆x, so that it becomes causally disconnected M˙BH (cid:20) (MBHβ1/2)(cid:21) from the ambient medium before encountering the sink M˙ 4.4 −1 ⊥ rveogluiomne. aFvoerraegaechofrunF, w=evc/avlFcuflaotreththeecemllassesi-twheeirgh5teodr M˙BH =min(1 , MBH(cid:20)1+ (MBHβ1/2)(cid:21) ) .(42) M 6 ∆x from the sink particle. Recall that gas is removed The perpendicular rate, however, potentially requires fromfluxtubesinsidethesinkregion. Theresultinglow- knowledge of both and β individually. Figures 9 density flux tubes are interchange unstable and will rise MBH and 10 plot these fits for the parallel and perpendicular awayfromtheaccretor. Sinceweareinterestedinverify- orientations. These fits are able to reproduce the simu- ing that the accreting gas is causally disconnected from lation data to within a factor of three. the ambientmedium, we include only accretinggas (i.e., For typical molecular cloud values of 5 and gas from cells where v x < 0, where x is the position M ∼ β 0.04 (Crutcher 1999), which corresponds to = vector from the sink par·ticle’s center) in calculating the ∼ MA 0.71, = 5.15, and β1/2 = 1.03, Equation averageof MF. This averagedvalue is given in Table 1. (33) gMiveBsHan accretion ratMe oBfH4.2 10−3 M˙ . For the For all our runs, we confirm that we have captured the B × fiducial parameters given in Equation (4), this corre- transition. With the simulation data, we can now determine βch sponds to 4.3×10−9 M⊙ yr−1, which is not that differ- andninourproposedinterpolationformulas(Equations entfromthehydrodynamicalprediction(Equation10)of 27 and 30) from Section 2.2 . Recall that these formu- M˙ =6.9 10−9 M yr−1. However,for smaller Mach BH ⊙ ×

Bondi-Hoyle Accretion in an Isothermal Magnetized Plasma

Similar Bondi-Hoyle Accretion in an Isothermal Magnetized Plasma