Feedback Effects on Low-Mass Star Formation

1.9 MB·English

by  Charles E. Hansen

The Astrophysical Journal, Accepted
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
FEEDBACK EFFECTS ON LOW-MASS STAR FORMATION
Charles E. Hansen 1, Richard I. Klein 1,2, Christopher F. McKee 1,3, and Robert T. Fisher 4
The Astrophysical Journal, Accepted
ABSTRACT
Protostellar feedback, both radiation and bipolar outflows, dramatically affects the fragmentation
and mass accretion from star-forming cores. We use ORION, an adaptive mesh refinement (AMR)
gravito-radiation-hydrodynamics code, to simulate low-mass star formation in a turbulent molecular
cloud in the presence of protostellar feedback. We present results of the first simulations of a star-
2 formingclusterthatincludebothradiativetransferandprotostellaroutflows. Werunfoursimulations
1 to isolate the individual effects of radiation feedback and outflow feedback as well as the combination
0 of the two. We find that outflows reduce protostellar masses and accretion rates each by a factor of
2 threeandthereforereduceprotostellarluminositiesbyanorderofmagnitude. Thismeansthat,while
radiation feedback suppresses fragmentation, outflows render protostellar radiation largely irrelevant
n
forlow-massstarformationaboveamassscaleof0.05M . Wefindinitialfragmentationofourcloud
a (cid:12)
J athalftheglobalJeanslength,around0.1pc. Withinsufficientprotostellarradiationtostopit,these
0.1 pc cores fragment repeatedly, forming typically 10 stars each. The accretion rate in these stars
3
scales with mass as predicted from core accretion models that include both thermal and turbulent
1
motions; the accretion rate does not appear to be consistent with either competitive accretion or
] accretion from an isothermal sphere. We find that protostellar outflows do not significantly affect
R the overall cloud dynamics, in the absence of magnetic fields, due to their small opening angles and
S poor coupling to the dense gas. The outflows reduce the mass from the cores by 2/3, giving a core
h. to star efficiency, (cid:15)core (cid:39) 1/3. The simulations are also able to reproduce many observation of local
star-formingregions. Oursimulationwithradiationandoutflowsreproducestheobservedprotostellar
p
luminosity function. All of the simulations can reproduce observed core mass functions, though we
-
o find they are sensitive to telescope resolution. We also reproduce the two-point correlation function
r of these observed cores. Lastly, we reproduce IMF itself, including the low-mass end, when outflows
st are included.
a
[
1. INTRODUCTION number of stars with some efficiency, 0.2 < (cid:15) < 1.0,
core
1 the IMF can be recreated. The actual conversion from
The origin of the stellar initial mass function (IMF)
v observed core masses to stellar masses may not be so
is one of the most fundamental problems of star forma-
1 simple due to cores blending in projection (Kainulainen
tion. The IMF can be described by single power law for
5 et al. 2009b; Michel et al. 2011), small cores that dis-
stars above 0.5 M (Salpeter 1955), and a broken power
7 (cid:12) perse before making stars (Myers 2009; Padoan & Nord-
2 law (Kroupa 2002) for stars below this mass. Alterna-
lund 2011) and cores accreting mass over time (Padoan
. tively, it can be described as a log-normal distribution
1 &Nordlund2011). Nevertheless, theCMFlikelyplaysa
with characteristic mass m = 0.2M that joins with
0 c (cid:12) strong role in creating the IMF
theSalpeterpowerlawforstarsabove1.0M (Chabrier
2 (cid:12) The observed CMF provides support to core accre-
2005). Any theory of the IMF must explain both the
1 tion theories of star formation (Shu 1977; McKee & Tan
functional form and the characteristic mass. A tantaliz-
: 2003), which start with a bound core and produce a sin-
v ing observational clue to the functional form lies in dust
gle stellar system. Simulations of turbulence find the
i observations in star-forming regions (Motte et al. 1998;
X functional form of the core mass function (log-normal
Testi&Sargent1998;Johnstoneetal.2000,2001;Motte
r et al. 2001; Beuther & Schilke 2004; Stanke et al. 2006; plus power law) is the expected outcome of supersonic
a turbulence (Padoan & Nordlund 2002; Padoan et al.
Alvesetal.2007;Enochetal.2008;Sadavoyetal.2010).
2007). Analytic predictions of a turbulent density field
These dust maps find many high density concentrations
with self-gravity can also reproduce this functional form
thatareconsistentwithprestellarandprotostellarcores.
(Hennebelle & Chabrier 2008, 2009). The characteristic
When the mass of these cores is calculated, the core
coremassisthentheJeansmassatsomecriticaldensity
mass function (CMF) has the same functional form as
and temperature. However, choosing the correct den-
the IMF, but with a higher characteristic mass, ranging
sity and temperature is problematic. In purely isother-
from0.2M to1M . Ifeachcoreisconvertedtoasmall
(cid:12) (cid:12) mal turbulence, there is no characteristic Jeans mass.
As objects collapse and the density increases, the Jeans
1Astronomy Department, University of California, Berkeley,
mass decreases. There is no transition where this de-
CA94720
2Lawrence Livermore National Laboratory, Livermore, CA creaseinJeansmasswillstopwithoutadditionalphysics.
94550 This means the core masses are either infinitely small or
3PhysicsDepartment, UniversityofCalifornia, Berkeley, CA functions of the global Jeans mass of the host molecu-
94720
4Physics Department, University of Massachusetts, Dart- lar cloud. Observations are consistent with a universal
mouth,MA IMF, however, even over a range of cloud Jeans masses
2
(Kroupa 2002; Chabrier 2003; Bastian et al. 2010). This
TABLE 1
means the characteristic core mass must be set by lo- Table of simulations
calphysics,whichisothermalturbulencecannotprovide.
Name ThermalPhysics Winds?
Star-forming regions are approximately isothermal be-
B Barotropic No
cause the thermal time scales are much shorter than the
BW Barotropic Yes
dynamical time scales, but there are ways to break this
R Radiation No
isothermality.
RW Radiation Yes
One approach is to focus on the coupling between
gasanddustinstar-formingenvironments(Larson2005; alytical estimates of mass loss from winds can fully ex-
Elmegreen et al. 2008). At low densities, gas-dust cou- plain the range of mass loss expected from observations
pling is poor and the gas is theoretically slightly sub- 0.2 < (cid:15) < 1.0 depending on the details of the cores
core
isothermal (temperature decreases with increasing den- andtheoutflows(Matzner&McKee2000). Theoutflows
sity). At higher densities, gas and dust are well coupled travelbeyondtheirstarsoforiginanddepositenergyinto
and the gas is theoretically slightly super-isothermal. parsec-scale turbulent motions. Evidence suggests that
This yields a critical density and temperature at the molecular cloud turbulence appears on the scale of the
transitionthatcanbeconvertedintoaJeansmass. This entire cloud (Ossenkopf & Mac Low 2002; Brunt et al.
criticaldensity,ρ∼10−19 gcm−3,islowerthantheden- 2009),soismostlikelydrivenbysourcesotherthanpro-
sities of large star-forming regions like Orion, however, tostellar outflows. Nonetheless, the dynamics on parsec
and unlikely to explain the characteristic core mass in scales can be significantly altered by outflows (Norman
these regions. &Silk1980;McKee1989;Li&Nakamura2006;Banerjee
One critical mass is the point when dust becomes et al. 2007; Nakamura & Li 2007; Swift & Welch 2008;
opaque to its own thermal radiation (Low & Lynden- Carroll et al. 2010; Arce et al. 2010; Wang et al. 2010).
Bell 1976). At that density, the gas will heat up and In this paper, we investigate the fragmentation of a
raise the Jeans mass, creating a minimum Jeans mass of parsec-scale molecular cloud into cores and then into
fragmentation. A barotropic simplification of this effect stars. This requires refinement to capture the fragmen-
sets the mass in many simulations (e.g. Bate & Bon- tation and radiative transfer to fragment at the correct
nell 2005; Bonnell et al. 2006; Offner et al. 2008; Bate mass scale, similar to Offner et al. (2009b). This also
2009a; Hennebelle et al. 2011). The density of this tran- requires simulation of protostellar outflows, to capture
sition is extremely high (∼ 10−13 g/cm3) (Masunaga the CMF to IMF efficiency, similar to the work of Cun-
et al. 1998) and the resulting Jeans mass (∼ 0.004M ) ningham et al. (2011) for high-mass star formation. The
(cid:12)
is much lower than the characteristic core mass (Low & goal of this paper is to explain both the observed CMF
Lynden-Bell 1976; Whitworth et al. 2007). In addition, and the IMF while self-consistently finding (cid:15) . In or-
core
the barotropic approximation is inaccurate when com- der to do this, we have performed the first simulation of
pared to simulations that include dust radiation (Boss a star-forming cluster to include both radiative transfer
et al. 2000; Krumholz et al. 2007a; Offner et al. 2009b; and protostellar outflows.
Bate2009b;Price&Bate2009;Tomidaetal.2010). The We describe our numerical method and simulation
importanceofdustradiationcanbeseeninBate(2009b) setup in §2. In §3, we report the results of our sim-
andPrice&Bate(2009),whofoundthattheinclusionof ulations. We discuss the implications of our results on
radiationsignificantlysuppressestheformationofbrown starformationtheoryandcomparetoobservationsin§4.
dwarfsdespitethenearabsenceofprotostellarluminosity We summarize our conclusions in §5.
in the simulations.
The most powerful break from isothermality comes 2. SIMULATIONS
from protostellar radiation. Massive protostars are ca-
Weperformfourprimarysimulationswithnearlyiden-
pable of heating an entire cloud (Krumholz et al. 2007a;
tical initial conditions but different physics as controlled
Cunningham et al. 2011; Myers et al. 2011; Krumholz
numerical experiments in order to isolate the impor-
et al. 2011). Low-mass stars do not have the same long
tance of feedback effects. These simulations all include
range influence, but simulations show they can still dra-
hydrodynamics, gravity and basic sink particle physics
matically reduce fragmentation in the disk and recover
(Krumholz et al. 2004), but may also include radiation
a 1 M characteristic core mass (Offner et al. 2009b;
(cid:12) (Krumholz et al. 2007a) and/or sink particle outflows
Krumholz et al. 2011). Protostellar radiation does not
(Cunningham et al. 2011). The simulations are shown
create a unique critical density, but it does weaken the
in Table 1. Barotropic simulations are labeled with a
density dependence of the effective Jeans mass (Bate
B, simulations with radiation are labeled with R, and
2009b).
simulation with protostellar winds are labeled with W.
Given a core mass function, there is still the question
Simulation labels can contain multiple letters.
of CMF to IMF efficiency, (cid:15) . The primary mecha-
core
nism for reducing the core mass is protostellar outflows. 2.1. Initial Conditions
Stars of all masses show bipolar outflows during their
All simulations have the same initial conditions, also
formation (Richer et al. 2000; Shepherd 2003) and are
used in Offner et al. (2009b). The initial gas tempera-
recreated in MHD simulations with sufficient resolution
ture is T = 10 K, the box length is L = 0.65 pc and
(Ciardi & Hennebelle 2010; Tomida et al. 2010). These g
outflows remove mass that would otherwise accrete onto the average density is ρ¯ = 4.46×10−20 g cm−3, corre-
stars, thereby reducing the final mass (Matzner & Mc- sponding to nH =1.91×104 cm−3. The total box mass
Kee 2000; Arce & Sargent 2006; Wang et al. 2010). An- is 185 M(cid:12). For radiative simulations, the radiation tem-
perature, T is initialized to 10 K. The radiation energy
r
3
density is thus E =aT4 =7.56×10−11 erg cm−3. dynamics including self-gravity, radiative transfer, pro-
r
Toobtaintheturbulentinitialconditions,webeginour tostellar outflows, and radiating star particles, all on an
simulations without self-gravity and apply velocity per- adaptive grid. Every cell in the grid has four conserved
turbationstoaninitiallyconstantdensityfieldusingthe quantities: mass density, ρ, vector momentum density,
method described in Mac Low (1999). These perturba- ρv, gas energy density, ρe, and radiation energy density,
tions correspond to a Gaussian random field with a flat E. These conserved quantities can be used to calculate
powerspectrumintherange1≤k ≤2. Theapplication derivedquantitiessuchasvelocity,v,andpressure,P. In
of these perturbations continues for three cloud cross- addition to the gas quantities, we evolve point-mass star
ing times and then stops. At this point the turbulence particles, each with a position x , mass M , momentum
i i
follows a Burgers power spectrum, P(k) ∝ k−2, charac- p , angular momentum, j and luminosity L . The sub-
i i i
teristic of supersonic hydrodynamic turbulence. The 3D script i refers to the star particle number. The particle
turbulent Mach number is M = 6.6, which gives a 3D methodisexplainedinKrumholzetal.(2004)(hereafter
rms velocity dispersion, σ =1.2 km/s. With this Mach KKM04),withtheadditionofradiation(Krumholzetal.
v
number the cloud is approximately virialized: 2007a) and outflows (Cunningham et al. 2011). The full
set of evolution equations for gas and particles is
5σ2
α = (cid:39)1. (1)
vir GM/R
∂ρ+∇·(ρv)+(cid:88)[M˙ W(r )−M˙ W (r )ξ(θ )]=0,
This is slightly above the linewidth-size relation ∂t KKM04 i w,i w i i
σ (cid:39) 0.7(R/1 pc)1/2 km s−1 (Solomon et al. 1987; i
(5)
Heyer & Brunt 2004), and is equivalent to σ =
1.2(R/1 pc)1/2 km s−1,whichiswellwithintheobserved
range (e.g. Falgarone et al. 2009) ∂ρv +∇·(ρvv)=−∇P −ρ∇φ−(cid:88)(p˙W(r )− (6)
After driving for three cloud crossing times, we then ∂t i
i
turn off driving, turn on gravity and follow the subse-
M˙ v W (r )ξ(θ )·ˆr ),
quentgravitationalcollapseforapproximatelyoneglobal w,i w,i w i i i
free fall time:
(cid:114)
3π ∂(ρe)
tff = 32Gρ¯ =0.315 Myr, (2) ∂t +∇·[(ρe+P)v]=ρv∇φ−κRρ(4πB−cE)−
where ρ¯is the mean density of the box. The simulations (cid:18) ρ (cid:19)2
Λ(T )− (7)
withradiationbecomeprohibitivelycomputationallyex- µm g
H
pensive at late times and are stopped at t = 0.83 t
ff (cid:88)
withatotalstellarmass of30M forsimulationR.The [ε˙ W(r )−
(cid:12) KKM04 i
barotropicsimulationsarecontinuedtot=1.05tff before i
they are stopped. At this time the total stellar mass in k T K
simulation B is 50 M(cid:12) compared to the total simulation M˙w,iWw(ri)ξ(θi)µ(Bγ−w1)],
mass of 185 M . There is still gas bound to protostars
(cid:12)
totaling 11 M when the simulations end. Our stellar
(cid:12) (cid:18) (cid:19)
mass estimates may therefore be too low by 20%. ∂ cλ
E−∇· ∇E =κ ρ(4πB−cE)+ (8)
Given our temperature of 10 K, the Jeans length at ρ¯ ∂t κ ρ P
R
is
λJ =(cid:18)πGcρ¯2s(cid:19)1/2 =0.20 pc, (3) (cid:18)µmρH(cid:19)2Λ(Tg)+(cid:88)LiW(ri),
i
and the Jeans mass is ∇2φ=−4πG[ρ+(cid:88)M δ(r )], (9)
i i
4π (cid:18)λ (cid:19)3 i
M = J ρ¯=2.7 M . (4)
J 3 2 (cid:12) 1
M˙ = M˙ , (10)
i 1+f KKM04
The turbulent Jeans mass, at density ρ = M2ρ¯, is 0.4 w
MT(cid:12)h.e calculations have a 2563 base grid with 4 levels of M˙w,i =fwM˙i = 1+fwf M˙KKM04, (11)
refinement by factors of 2, giving an effective resolution w
of 40963. This resolution corresponds to ∆x4 = 32 AU p˙i =p˙KKM04, (12)
at the finest refinement level.
r =x−x , (13)
i i
2.2. Evolution Equations
θ =acos(ˆr ·ˆj). (14)
We use the parallel adaptive mesh refinement code i i i
ORION for our simulations. The numerical method is The quantities entering these equations are defined
nearly identical to what we have used in previous pa- below. Equations (5) and (6) are the fluid equations
pers (Krumholz et al. 2007a; Offner et al. 2009b; Cun- for mass and momentum, modified to include particles.
ningham et al. 2011; Krumholz et al. 2011; Myers et al. Equations (7) and (8) are the energy equations for gas
2011). ORION solves the equations of compressible gas and radiation respectively. The Poisson equation for the
4
gravitational potential, φ is given by equation (9). The overthe4-cellaccretionzoneoftheparticleandW rep-
w
particle evolution is given by equations (10), (11) and resentstheoutflowwindowfunction. Outflowsareimple-
(12). We use periodic boundary conditions for all gas mentedasinCunninghametal.(2011),andsummarized
and particle quantities. here.
For the radiative runs, we adopt Marshak boundary Our outflow model is specified by the dimensionless
conditions for the radiation field. This allows radiation parameter f , which sets the mass flux of outflow as
w
to escape from the box as it would from a molecular a fraction of the accretion onto a star, and v , the
w,i
cloud. The equation of state for the gas is given by wind launch speed. The mass fraction in our simula-
tions is f = 0.3. The wind speed is set by the Keple-
ρk T (cid:18) v2(cid:19) w (cid:112)
P = B g =(γ−1)ρ e− , (15) rian speed at the surface of the star, vk,i = GMi/r∗,i
µmH 2 where r∗,i is the protostellar radius, but is is capped
at 60 km s−1 for computational speed. Specifically,
where µ=2.33 is the mean molecular weight for molec-
v =min(v ,60 km s−1). The velocity cap has a sim-
ular gas of Solar composition and γ is the ratio of spe- w,i k,i
ilar effect to the choice of Cunningham et al. (2011) to
cific heats. Since most of the simulation domain is too
use v = v /3 for the most massive stars in the cal-
coldtorotationallyexcitemolecularhydrogen, weadopt w,i k,i
γ = 5/3, representing a monatomic ideal gas. The term culation. Wind gas is injected at temperature Tw =104
κ ρ(4πB−cE) in equations (7) and (8) represents en- K.
P
ergy exchanged between the radiation field and the dust The wind is injected over a window function Ww,
in our gas, with B = caTg4/4π representing the Planck 1 (cid:26)r−2 if 4∆x<r ≤8∆x
emission function integrated over all frequencies. The W = . (17)
w C 0 otherwise
opacities κ and κ are Planck and Rosseland means 1
P R
given by the dust opacities from the iron normal, com-
which represents a shell just outside of the accretion re-
posite aggregates dust model of Semenov et al. (2003).
gion. The normalization constant C is computed nu-
1
We assume that the gas and the dust are thermally cou-
merically to avoid numerical aliasing effects that occur
pled. When the gas temperature exceeds the dust de-
from injecting a spherical wind into a Cartesian grid.
struction temperature, the energy exchange term goes
Theexactangulardistributionofthewindisdescribed
to zero and the gas and radiation unrealistically decou-
in the function ξ. The functional form is taken from
ple. Toaddresscoolingfromgasabovethedustdestruc-
Matzner & McKee (1999),
tion temperature, we use the line cooling function Λ(T )
g
from Cunningham et al. (2006). This removes energy (cid:20) (cid:18) 2 (cid:19) (cid:21)−1
from the gas and adds that energy to the radiation field ξ(θ)= ln (sin2θ+θ2) , (18)
θ 0
(see Cunningham et al. (2011) for further details). The 0
radiation flux limiter is given by λ = 1 (cid:0)cothR− 1(cid:1), where θ is a flattening parameter that sets the opening
R R 0
whereR=|∇E|/κ ρE (Levermore&Pomraning1981). angleofthewind. Weusethefiducialvalueofθ =0.01.
R 0
It should be noted that we have excluded the radiation Equation(18)isaveragedoverthesolidanglesubtended
pressure and radiation enthalpy advection terms from by a grid cell at the outer radius of W . This averaging
w
equations (6), (7) and (8) that appear in the analogous is particularly important near θ = 0. In addition, ξ is
equationinKrumholzetal.(2007a). Thisapproximation set to zero for θ ∼ π/2, so that winds are not injected
is justified in the formation of low-mass stars, as shown directly into the plane of any equatorial disks.
in Offner et al. (2009b). We update the luminosity of each star, L , using the
i
Whenradiationisneglected,theenergyexchangeterm protostellar evolution model described in Offner et al.
fromequation(7)disappears,andweclosethesystemof (2009b). Inthismodel,75%oftheaccretionenergyisra-
equations with a barotropic equation of state for the gas diatedawaywhile25%isnominallyusedtopowerwinds.
pressure: The energy of winds in our simulations is already deter-
(cid:34) (cid:18) ρ (cid:19)γ−1(cid:35) minedbyfw andvw,i,whichistypically15%oftheaccre-
P =ρc2 1+ , (16) tion energy. The remaining 10% of the accretion energy
s0 ρ
c is effectively lost. When winds are not present, we still
only radiate 75% of the accretion energy for consistency
(cid:112)
wherecs0 = kBT0/µmH istheisothermalsoundspeed across simulations.
at temperature T0 = 10 K, γ = 5/3 and ρc is the crit- The evolution equations can be described as the fluid
ical density. The critical density determines the transi- and radiation equations from Offner et al. (2009b) com-
tion from isothermal to adiabatic regimes and we adopt binedwiththeparticleequationsandlinecoolingofCun-
ρc = 2×10−13 g cm−3 to agree with the collapse so- ningham et al. (2011), but there is one important mod-
lution from Masunaga et al. (1998) prior to H2 dissoci- ification. In the KKM04 sink particle methodology, all
ation. Simulations that use the barotropic equation of particleswithoverlappingaccretionzonesaremergedto-
state achieve maximum densities of 5×10−15 g cm−3, gether. Thisgivesaneffectivemergerradiusof8cells,or
with an effective temperature at ρmax of 10.8 K. Most of 256AUatagridresolutionof32AU.Tolimitthiseffect,
thegasinanygivensimulationsiseffectivelyisothermal. wechangedthemergerradiusto4cells,representingthe
The particle quantities M˙ , p˙ and point when a particle is in the accretion zone of another
KKM04 KKM04
ε˙ in Equations (5 - 7) represent the sink parti- particle. This gives an effective merger radius of 128 AU
KKM04
cle accretion rates of mass, momentum and energy from at our resolution. Even with this improvement, our par-
the surrounding gas in the absence of winds as given by ticle algorithm will unrealistically merge stars that pass
KKM04. The function W represents a window function within 128 AU. To address this, we have implemented a
5
mass limit of m = 0.05 M , above which stars do
merge (cid:12)
notmerge. Thislimitischosentocorrespondtothemass
of second collapse in the formation of a star with final
mass1M (Masunagaetal.1998;Masunaga&Inutsuka
(cid:12)
2000). Second collapse occurs when the protostar’s core
temperature becomes high enough to dissociate molec-
ular hydrogen. Before second collapse, protostars are
extended balls of gas with radii of a few AU and have a
muchhighercollisionalcross-sectionthanmainsequence
stars. Aftersecondcollapse,stellarmergersshouldbeex-
tremelyrareandwedonotallowthem. Thisapproachis
alsousedinMyersetal.(2011). Themassofsecondcol-
lapse depends on the accretion history of the protostar
and is necessarily less than 0.05 M for brown dwarfs
(cid:12)
(Stamatellos&Whitworth2009;Bate2011). Theeffects
of numerical merger suppression is explored in §4.9.
ORION utilizes a second order Godunov scheme to
solve the equations of compressible gas dynamics (Tru-
elove et al. 1998; Klein 1999). These are equations (5)-
(7), excluding terms from stars and radiation. The Pois-
son equation (9) is solved using a multi-grid iteration
scheme (Truelove et al. 1998; Klein 1999; Fisher 2002).
Theflux-limiteddiffusionradiationequation(8)andthe
radiation terms in equation (7) are solved using the con-
servative update scheme from Krumholz et al. (2007b)
modified to include the pseudo-transient continuation of
Shestakov & Offner (2008).
We use the Truelove criterion (Truelove et al. 1997) to
determinetheadditionofnewAMRgridssothatthegas
density in the calculation always satisfies
J2πc2
ρ< s , (19)
G(∆x )2
l
where ∆x is the cell size on level l. We adopt a Jeans
l
numberof0.125. Inthesimulationswithradiativetrans-
fer, it is necessary to resolve the spatial gradients in the
radiationfield. Areasofhighradiationgradientsarenear
accreting stars, which tend to already be refined under
the Truelove criterion. This is not always true for more
evolved stars, which have higher luminosities and have
accretedthedensegasthatwouldtriggerrefinement. We
find that the radiation gradients are adequately resolved Fig. 1.—ColumndensityoftheentiresimulationdomainforBW
(left)andB(right)attimes0,0.2,0.4,0.6,0.8and1.0t fromtop
by refining whenever |∇E|∇x /E >0.25. ff
l to bottom. Star particles are marked with white circles. There is
very little difference on the domain scale with and without winds
3. RESULTS for the barotropic simulations. The color bar is gcm−2 and the
entiredomainis0.67pcacross.
3.1. Large Scale Evolution
The evolution of the barotropic simulations with and
without winds is depicted in Figure 1. Figures 2 and 3 in the core and an additional group of low-mass stars
depict the evolution of the radiative simulations without forms, totaling ∼ 10 stars per core. These cores with
and with winds, respectively. In all simulations, for multiple stars generally resemble observed high-stellar-
t(cid:46)0.4t ,therearecloud-scalefilamentsthatslowlycon- density cores (Teixeira et al. 2007; Chen & Arce 2010).
ff
tract, allowing 3 turbulent cores of width ∼ 20,000 AU There are an additional 20 stars in cores that are still
and density ∼10−19 g cm−3 to form. This length is half forming at the end of the simulation, giving a global to-
the Jeans length at the average cloud density. The first tal of 80 stars. Three of the cores coalesce by the end of
coretoformisatthecenterofeachpanelinFigures1-3, the simulation to form a single group of 30 stars.
the second core is left of the central core, near the left The evolution of the 3D rms velocity dispersion, σ , is
v
edge,andthethirdcoreisneartherightedge,attheend shown in Figure 4. The global turbulence decays un-
ofasimulation-widefilament. Atthispoint,thecoresbe- til star formation ramps up at t ∼ 0.5t . There are
ff
gintofragment,whilenewcoresform,eventuallyforming two main mechanisms for star formation to increase σ .
v
6 fully developed cores. These cores each have a cen- First, as stars accrete mass and deepen their gravita-
tral stellar system. In simulation R, this central system tionalpotential,thesurroundinggascanconvertgravita-
is all that forms in each core. In all other simulations, tionalenergyintokineticenergyasitfallsintothestars.
the central system represents ∼ 75% of the stellar mass This is shown in the gradual increase in Mach number
6
Fig. 2.—Columndensity(left)anddensityweightedtemperature Fig. 3.— Same as Figure 2, but for simulation RW. The high
(right) for simulation R at times 0, 0.2, 0.4, 0.6, and 0.8 t from temperature regions are the paths of outflows. It only takes a
ff
toptobottom. Starparticlesaremarkedwithwhitecircles. small amount of gas at 104 K to move the average temperature
above30Kthroughthatlineofsight.
for t > 0.6t in the B simulation. This effect is strong
ff
injected by the winds is greater than the characteris-
enough to return σ to near its original virialized value
v
tic cloud momentum. At this point, the total amount
by itself. In rare cases, a many-body close encounter be-
of turbulent momentum that has been dissipated (in-
tween stars will eject some gas at high velocities. There
cluding dissipation of wind momentum) is roughly the
is not much momentum injected this way and the en-
total amount of momentum that has been injected by
ergy quickly dissipates, but it causes spikes in σ for the
v
the winds. By the end of the BW and RW simulations,
barotropicsimulations,whichhavemoresmall-scalefrag-
the total wind momentum injected into the cloud is over
mentation and therefore more many-body close encoun-
twice the characteristic cloud momentum. The kinetic
ters. The second mechanism occurs when protostellar
energyinjectedfromthewindsdissipatesovertime,sug-
winds are included. Somemass accreted ontostars isdi-
gesting a steady-state solution where the velocity dis-
rectly injected around the stars at high velocities. This
persion of the cloud is constant with time as the winds
causesthesmoothincreaseinσ forsimulationsBWand
v
replenish energy as quickly as it can dissipate.
RW as well as spikes from events with particularly high
accretion rates that lead to bursts in wind momentum.
3.2. Evolution of the Protostellar Population
ThetotalmomentuminjectedbywindsformodelBW
is shown in Figure 5. For comparison, a characteris- The total mass in stars and the total number of stars
tic value of the magnitude of the momentum associated as functions of time is shown are Figures 6 and 7. The
withinternalmotionsinthecloud,M σ isalsoplot- realizationoftheinitialturbulenceisslightlydifferentbe-
cloud v
ted. For t > 0.8t , the total momentum that has been tween the radiative and barotropic simulations, so star
ff
7
however, because protostellar luminosity inhibits frag-
mentation and the winds reduce that luminosity. The
threesimulationsotherthanRshowadramaticincrease
in the number of stars at 0.6 t <t<0.8 t .
ff ff
Fig. 4.—TimeevolutionofglobalrmsMachnumberforsimula-
tionswithandwithoutwindsandwithandwithoutradiation. The
turbulentenergyinallsimulationsdecaysforhalfaglobalfreefall
time, at which point gravitational potential energy from stars is
convertedintokineticenergy,whichraisesthermsvelocity. When
winds are included, they contribute over twice as much energy as
gravityitself.
Fig. 6.— Total mass in stars as a function of time for the four
simulations. Themassoftheentiresimulationdomainis180M(cid:12).
Fig. 5.—Totalmomentumthathasbeeninjectedbywindsover
timeforthebarotropicsimulations. Forcomparison,thetotalmass
ofthesimulationmultipliedbythevelocitydispersionisalsoplot-
ted. Thetotalwindmomentumintegratesallinjectedmomentum
over time, even from winds that have decayed. As a result, the
amountofinjectedmomentumiseventuallyhigherthantheactual
momentumofthecloud.
formation begins at different times. This is the result
of a slightly higher maximum density due to turbulence
anditisunrelatedtothechoiceofradiativeorbarotropic
thermodynamics. Thefirststarsformatt=0.2t forthe
ff
radiative simulations and t = 0.35t for the barotropic
ff
simulations. The turbulence overlaps enough between
theradiativeandbarotropiccases,however,thatatlater
times the total mass in stars at any given time is similar
for the two cases when winds are not included. Winds
reduce the mass in stars by about a factor of 3 in both
the radiative and barotropic cases. The number of stars
Fig. 7.—Numberofstarsasafunctionoftimeforthefoursim-
does not change between BW and B, implying winds do
ulations. Inthebarotropiccase,thenumberofstarsisunaffected
not cause or suppress fragmentation by themselves. The by winds. In the radiative case, radiation suppresses the number
number of stars in RW is significantly greater than in R, ofstarsunlesswindsarepresent.
8
The evolution of median stellar mass is shown in Fig-
ure8. Thisisaroughproxyforthecharacteristicmassof
the protostellar mass function. Note that it will always
be lower than the median mass of the IMF, because not
all of the stars have finished accreting. Both BW and
RW maintain a median around 0.05 M (similar to that
(cid:12)
of the protostellar mass functions in McKee & Offner
2010) throughout the simulation. The median does in-
creaseforsimulationBWaroundt>t astheformation
ff
rate of new stars decreases. The median of B fluctuates
more, but is around 0.2 M . Lastly, R maintains a me-
(cid:12)
dian around 0.5M . This general behavior should be
(cid:12)
expected. Themedianmassislowestwhenwindsarein-
cludedandhighestwhenradiationisallowedtosuppress
fragmentation. The case with both winds and radiation
endsupsimilartoBWbecausewindsreduceprotostellar
luminosities.
Fig. 9.—Totalstellarluminosityversustimeforsimulationswith
radiation,bothwithandwithoutwinds. Windsdramaticallylower
theluminosity.
As was seen in the plot of total luminosity, the mean
and median values of protostellar luminosity are much
lower when winds are included. The disparity in aver-
age luminosity is even greater than the disparity in to-
tal luminosity because there are fewer stars when winds
are excluded. The average luminosity in simulation R is
heavily influenced by a single 6.6 M star that accounts
(cid:12)
for over half the total luminosity in the simulation. Un-
like the low-mass stars, most of this luminosity is pow-
eredbynuclearfusionratherthanaccretion. Protostellar
luminosities will be discussed further in section 4.3
Fig. 8.— Median mass of stars as a function of time for the
foursimulations. Thetwocaseswithwindsmaintainlowmedians
throughoutthesimulations. Thecasewithradiationwithoutwinds
(caseR)isabletosuppressfragmentationandnewstarformation
largelystopsastheoriginalstarsaccretemass.
Thegloballuminosityevolutionfortheradiativesimu-
lationsisplottedinFigure9. Thewindsreducethetotal
luminositybyafactorofupto10atanygiventime. This
is expected since the global accretion luminosity is
L= 3 (cid:88)GM(cid:63)iM˙(cid:63)i ∼ GM(cid:63),totM˙(cid:63),tot; (20)
4 R (cid:104)R (cid:105)
(cid:63)i ∗
i
given that the total mass in stars M is reduced by
(cid:63),tot
a factor of 3 when winds are included, and the total
accretion rate of stars M˙ is also reduced by 3, the Fig. 10.— Mean and median stellar luminosity versus time for
(cid:63),tot
total luminosity is therefore reduced by a factor of 9 as- simulationswithradiation. Thetoppanelisthesimulationwithout
windsandthebottompanelisthesimulationwithwinds.
suming the characteristic stellar radius, (cid:104)R (cid:105), does not
(cid:63)
change. Main sequence stars typically have a positive
correlation between mass and radius, suggesting the fac-
3.3. Thermal Evolution
tor of 9 is an upper limit. However, at any given point
in time, many stars in our simulations are in the de- All simulations start at a background temperature of
generate regime where mass and radius are negatively 10 K and are bathed in 10 K radiation. Stellar radiation
correlated (Chabrier et al. 2009), which counteracts the and mechanical energy from protostellar outflows, can
positivecorrelationfromthehighermassstarsandkeeps raise this temperature. We have identified gas heated
the total luminosity ratio near 9. above 12 K as thermally affected by stellar feedback.
The average stellar luminosity is shown in Figure 10. Turbulent dissipation alone heats almost no gas above
9
12K,leavingstellarfeedbackastheonlyexplanationfor
this heating. The total mass of this gas is shown in Fig-
ure 11. The simulation with winds has significantly less
Fig. 11.—Totalgasmassheatedabove12Kversustime, com-
paredtothebackgroundvalueof10K.Thetotalmassinstarsis
alsoplottedforreference.
Fig. 12.— Phase plot showing total gas mass as a function of
temperature(y-axis)anddensity(x-axis)forradiativesimulations
heated gas than the simulation without winds. This is
with (left) and without (right) winds. Phase plots are taken at
duetothereducedluminositycausedbythewindsshown times of 0.25, 0.5 and 0.75 t from top to bottom. The high-
ff
in Figure 9. In each simulation, the mass in heated gas temperature, low-density gas on the left part of the wind phase
plotsisoutflowgas. Warm,high-densitygasisgasnearaluminous
roughly follows the mass in stars. When winds are in-
star.
cluded, the mass in stars drops by a factor of 3, and the
mass of heated gas also falls due to the reduced stellar
luminosity. Someofthelostradiativeheatingisreplaced mentally different from the isotropic, homogeneous hy-
by mechanical heating from outflows. This is evident in drodynamic turbulence it replaces. This result is also
the ratio of heated gas mass to stellar mass. The mean seen in Nakamura & Li (2007), who find that the late
value of this ratio is 0.9 without winds, but rises to 1.5 timeturbulentstatisticsdonotmatchexpectedisotropic
with them. Wind gas is injected at a temperature of 104 hydrodynamic results. One key difference is the energy
K,butcoolsquickly. Themassingaswithtemperatures from outflows is highly anisotropic. Outflow cavities are
above 1000 K is less than 2% of the mass in stars at any marked by long channels with high velocity shear be-
given time. tweenthefastoutflowgasandtheslowambientgas. This
Tofurtherexploretheheatedgas,temperature-density shear is detectable as solenoidal energy. There is some
phase plots with and without winds are shown in Figure compressiveenergyattheheadoftheoutflowcavity,but
12. The phase plots with and without winds are notably most of the surface area is the side walls of the cavity
different in two areas. First, the wind gas fills the high- and not the head. To measure the relative importance
temperature,low-densitydomain,whilethesamedomain of solenoidal and compressive motions, we use the ratio
is empty without winds. Second, high-density gas with (cid:104)|∇×v|2(cid:105)/(cid:104)(∇·v)2(cid:105), shown in Figure 13. In the case of
ρ>10−16g/cm3hasahighertemperaturerangewithout isotropicturbulence,thisratiois2,whichisalsotheratio
ofsolenoidaltocompressiveenergies(Elmegreen&Scalo
windsthanwithwindsbecausetheextrastellarluminos-
2004). Wind injection greatly increases the solenoidal
ity heats that gas. When winds are included, that dense
velocities, steadily increasing the solenoidal to compres-
gas is less common in addition to being colder; there is
sive ratio over the course of the simulation. Bursts in
morefragmentation,whichturnsdensegasintostars. In
outflows injection around t = 0.5t and t = 0.65t are
addition,someofthegasisalsoblownawaybythewinds ff ff
also visible in the solenoidal velocity. Anisotropic tur-
themselves.
bulence behaves differently than isotropic turbulence; in
4. DISCUSSION particular, it takes longer to decay since it decays on
the crossing time calculated from the smallest velocity
4.1. Supporting a Cloud with Outflows
dispersion (Hansen et al. 2011). It is difficult to mea-
The turbulent evolution shown in Figures 4 and 5 sure this increase in the decay time in our simulations,
roughly agrees with previous simulations of molecular however, becausethewindsdriveturbulenceoverawide
clouds with outflows, such as those in Nakamura & Li range of scales.
(2007) and Wang et al. (2010). Turbulent energy decays The other major difference between the initial turbu-
initially,onlytobereplacedbykineticenergyfromwinds lence and wind-driven turbulence is seen in the rms ve-
and from gravity. While these sources can increase the locity, σ , of gas with ρ > ρ¯. Even in isotropic,
dense
total kinetic energy of gas, the new turbulence is funda- homogenous, hydrodynamic turbulence, there is a neg-
10
a centrally concentrated cloud can achieve a quasi-static
balance between outflows and gravity.
4.2. Comparison to the IMF
Our simulations are marginally able to capture bina-
ries, so most star particles should represent a single stel-
larobjectinsteadofastellarsystem. Thetypicalbinary
separation for main sequence stars is 50 AU (Mathieu
1994), just slightly larger than our resolution of 32 AU.
We cannot form stars within 128 AU of another star
due to merging star particles, but stars that form be-
yond 128 AU can spiral in to 32 AU through interaction
with gas. At any given time, about 2/3 of our star par-
ticles are in stellar multiples. The multiple properties
are dynamic due to unstable high-order multiples and
we cannot compare to the observed system properties.
We can, however, compare to observed stellar proper-
ties. The mass functions of the four main simulations
Fig. 13.— Ratio of solenoidal to compressive velocities squared are shown in Figure 15 and compared to the stellar IMF
((cid:104)|∇×v|2(cid:105)/(cid:104)(∇·v)2(cid:105)) versus time with and without winds for in Chabrier (2005) as well as protostellar mass functions
the radiative simulations. For pure hydrodynamic, isotropic tur- from McKee & Offner (2010). The mass functions in
bulence,thisratioisaround2. Thisratiostaysnear2whenwinds
are excluded. When winds are included, the turbulence is much
moreanisotropic,leadingtohighersolenoidalfractions.
ative correlation between density and velocity, causing
σ <σ (Offner et al. 2009a). The winds themselves
dense v
are collimated, very low density gas and have difficulty
transmitting energy into high density gas. This means
that while σ is much greater with winds than without,
v
σ does not change much when winds are included.
dense
The evolution of σ compared to σ is shown in Fig-
dense v
ure 14.
Fig. 15.—Themassfunctionofallstarsineachsimulationare
shown in blue histograms. The stellar IMF from Chabrier (2005)
is plotted as the solid green line. The protostellar mass functions
fortheturbulentcoremodelandtheisothermalspheremodelfrom
McKee&Offner(2010)arethedashedlinesanddashdottedlines,
respectively. Top left: RW at t = 0.83t . Top right: R at t =
ff
0.83t . Bottom left: BW at t = 1.09t . Bottom right: B at
ff ff
t=1.03t
ff
Figure 15 are shown at the latest time available for each
simulation. The barotropic simulations could be evolved
to later times due to the computational expense of flux-
limited diffusion with many stellar sources. The mass
Fig. 14.—TimeevolutionofrmsMachnumberofdensegasand
functions in Figure 15 are not exactly comparable to an
allgaswithandwithoutwinds. WindssignificantlyraisetheMach
number of the light gas, but do not strongly influence the dense IMF because some of the stars are still accreting. The
gasturbulence. barotropic mass functions have evolved to a later time
and should be similar to the IMF. The radiative simu-
Because the dense gas is relatively unaffected by the lations are still actively accreting and are closer to the
outflows, if our cloud had been centrally concentrated protostellar mass functions in McKee & Offner (2010).
likethatofNakamura&Li(2007)orWangetal.(2010), Among the protostellar mass functions, the turbulent
the dense part of the cloud would have most likely col- core and competitive accretion (not shown in Figure 15)
lapsedonitselfevenwiththesupportofprotostellarout- models both roughly match RW. The isothermal sphere
flows. Magnetic fields may have an effect, as shown by protostellar mass function is the best fit to mass func-
Wang et al. (2010). The magnetic fields help transmit tionsofthesimulationswithoutwinds,butthisismerely
outflowenergytoamuchlargersolidangle, sothateven a function of both mass functions being too top-heavy.