Simulating star formation in Ophiuchus

1.5 MB·

by  O. Lomax

Mon.Not.R.Astron.Soc.000,1–??(2013) Printed31January2014 (MNLATEXstylefilev2.2)
Simulating star formation in Ophiuchus
O. Lomax(cid:63)1, A. P. Whitworth1, D. A. Hubber2,3, D. Stamatellos4 and S. Walch5,6
1School of Physics and Astronomy, Cardiff University, Cardiff CF24 3AA, UK
2University Observatory, Ludwig-Maximilians-University Munich, Scheinerstr.1, D-81679 Munich, Germany
3Excellence Cluster Universe, Boltzmannstr. 2, D-85748 Garching, Germany
4Jeremiah Horrocks Institute, University of Central Lancashire, Preston, Lancashire, PR1 2HE, UK
4 5Physikalisches Institut der Universita¨t zu K¨oln, Zu¨lpicher Strasse 77, D-50937 Cologne, Germany
1 6Max-Planck-Institut fu¨r Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching
0
2

n 31January2014
a
J

0 ABSTRACT
3 Wehavesimulatedstarformationinprestellarcores,usingSPHandinitialconditions
informedbyobservationsofthecoresinOphiuchus.Becausetheobservationsarelim-
]
itedtotwospatialdimensionsplusradialvelocity,wecannotinferinitialconditionsfor
R
the collapse of a particular core. However, with a minimum of assumptions (isotropic
S
turbulence with a power-law spectrum, a thermal mix of compressive and solenoidal
.
h modes,acriticalBonnor-Ebertdensityprofile)wecangenerateinitialconditionsthat
p match,inastatisticalsense,thedistributionsofmass,projectedsizeandaspectratio,
- thermalandnon-thermalone-dimensionalvelocitydispersion,observedinOphiuchus.
o
The time between core-core collisions in Ophiuchus is sufficiently long, that we can
r
t simulate single cores evolving is isolation, and therefore we are able to resolve masses
s
a well below the opacity limit. We generate an ensemble of 100 cores, and evolve them
[ withnoradiativefeedbackfromthestarsformed,thenwithcontinuousradiativefeed-

back, and finally with episodic radiative feedback. With no feedback the simulations
2
producetoomanybrowndwarfs,andwithcontinuousfeedbacktoofew.Withepisodic
v
radiative feedback,boththe peak oftheprotostellar massfunction (at 0.2M )and
7 ∼ (cid:12)
3 the ratio of H-burning stars to brown dwarfs are consistent with observations. The
2 mass of a star is not strongly related to the mass of the core in which it forms. Low-
.7 cmoarsessc(Mores (M(cid:38)co1rMe ∼)0u.1suMa(cid:12)lly) tfreangdmteonctoilnlatopsseevinertaolsoinbgjelectosbwjeitchts,awbhroearedams haisgshr-amnagses.
1 core (cid:12)
0
Key words: Stars:formation,stars:low-mass,stars:massfunction,stars:pre-Main-
4
Sequence, hydrodynamics.
1
:
v
i
X
r 1 INTRODUCTION In Turbulent Fragmentation, molecular clouds do not
a collapse on a global scale, but turbulent flows within the
The origin of the stellar initial mass function (IMF), and
cloud produce dense cores (e.g. Padoan & Nordlund 2002;
whyitappearstobeapproximatelyuniversal,remainmajor
Hennebelle&Chabrier2008,2009),roughlybetween0.01pc
problems in star formation. Two main scenarios have been
and 0.1pc across. Cores are usually well separated from
proposed to addresses these problems.
one another, and those that are gravitationally bound col-
InCompetitiveAccretion,partofamolecularcloudcol-
lapse to form stars (e.g. Andre, Ward-Thompson & Bar-
lapses, thereby providing a gas-rich environment in which
sony 1993, 2000). Star formation in cores has been sim-
stars form. The most massive stars preferentially accrete
ulated numerically by Bate (1998, 2000); Horton, Bate
the most gas; the least massive stars may be ejected from
& Bonnell (2001); Matsumoto & Hanawa (2003); Good-
thesystembyN-bodyinteractionsandceaseaccreting(e.g.
win & Whitworth (2004); Delgado-Donate, Clarke & Bate
Bonnell et al. 1997; Reipurth & Clarke 2001; Bonnell, Vine
(2004);Delgado-Donateetal.(2004);Goodwin,Whitworth
& Bate 2004; Bate & Bonnell 2005; Bate 2009a). This sce-
&Ward-Thompson(2004,2006);Walchetal.(2009,2010);
narionaturallyyieldsadistributionofstellarmassessimilar
Walch, Whitworth & Girichidis (2012).
in shape to the IMF (Kroupa 2001; Chabrier 2005), with a
peak close to the cloud’s original Jeans mass. Recent millimetre and submillimetre observations of
starformingregions(e.g.Motte,Andre&Neri1998;Testi&
Sargent1998;Johnstoneetal.2000;Motteetal.2001;John-
(cid:63) E-mail:[email protected] stoneetal.2001;Stankeetal.2006;Enochetal.2006;John-
c 2013RAS
(cid:13)
2 O. Lomax, A. P. Whitworth, D. A. Hubber, D. Stamatellos and S. Walch
stone&Bally2006;Nutter&Ward-Thompson2007;Alves, Ophiuchusaregroupedintoclumps1,labelledOph-A,Oph-
Lombardi&Lada2007;Enochetal.2008;Simpson,Nutter B1, Oph-B2, Oph-C, Oph-D, Oph-E and Oph-F.
& Ward-Thompson 2008; Rathborne et al. 2009; Ko¨nyves We base the initial conditions for our simulations on
etal.2010)haveshownthatthecoremassfunction(CMF), the observations of MAN98 and ABMP07. MAN98 present
issimilarinshapetotheIMF,butshiftedupwardsinmass continuum observations of cores with a completeness limit
byafactorof3to5.Fromthisitiscommonlyinferredthat of ∼ 0.1M . These data are complemented by observa-
(cid:12)
the shape of the IMF is directly inherited from the CMF, tionsofcorevelocitydispersionsinOphiuchusbyABMP07.
i.e. each core collapses and converts roughly the same frac- There are more recent measurements of core masses and
tionofitsmassintothesamenumberofprotostars.Holman sizes in Ophiuchus (e.g. Stanke et al. 2006; Simpson, Nut-
et al. (2013) have shown that this is only compatible with ter & Ward-Thompson 2008), but they contain fewer cores
the observed stellar multiplicity statistics if a typical core with an associated ABMP07 velocity dispersion than the
spawns 4 or 5 stars. MAN98 data. Further details of the MAN98 and ABMP07
The L1688 cloud (hereafter Ophiuchus) is a region observations are given in the next two sections.
wherestarsareformingincores.ObservationsofOphiuchus
by Motte, Andre & Neri (1998) and Andr´e et al. (2007)
2.1 MAN98
(hereafter MAN98 and ABMP07) show that the cores are
wellseparatedandunlikelytointeractdynamicallywithone MAN98 map dust emission at 1.3mm using the MPIfR
another before they collapse. In this paper we use SPH to bolometerarrayontheIRAM30mtelescope.Thebeamsize
simulate star formation in these cores, using initial condi- is 11(cid:48)(cid:48), corresponding to 1500AU at Ophiuchus, and they
tions informed as closely as possible by the MAN98 and aresensitivetohydrogencolumn-densitiesN>1022Hcm−2.
ABMP07observations.Alimitationofourapproachisthat Theirmapcoversanareaof∼480arcmin2 andincludesthe
cores tend to be embedded in larger scale clumps and fil- Oph-A,Oph-B1,Oph-B2,Oph-C,Oph-D,Oph-EandOph-
aments, and material that would accrete onto the cores as F clumps. MAN98 use a multi-scale wavelet analysis to ex-
they collapse is not included in the simulations. tract 61 cores, 36 of which are resolved. Two dimensional
Deriving reliable initial conditions for simulations of Gaussian distributions are fitted to the cores and the core
corecollapse isa difficultinverseproblem. Coremasses can dimensionsaregivenbytheFWHMofthemajorandminor
be estimated in a relatively straightforward way from ob- axes.Thecoreshavemassesbetween∼0.1M and∼3M ,
(cid:12) (cid:12)
servedsubmillimetrecontinuumfluxes.However,theintrin- and dimensions between ∼1000AU and ∼20,000AU.
sic three-dimensional sizes and shapes of cores can only be Thecoremassesandsizesusedinthispaperareslightly
constrainedbyfittingtheparametersofacrediblemodelto differentfromthosegiveninMAN98.Themassofaprestel-
observed data (see Lomax, Whitworth & Cartwright 2013, lar core is calculated using
hereafter LWC13). The detailed three-dimensional motions
S(λ)D2
within cores are even less well constrained by observations M = , (2.1)
of radial velocity profiles. A key element of this paper is to core κ(λ)B(λ,T)
generate an ensemble of cores which, in a statistical sense, whereS(λ)isthemonochromaticfluxfromthecoreatwave-
resemble those observed in Ophiuchus. lengthλ,Disthedistancetothecore,κ(λ)isthemassopac-
In Section 2 we review the data used to constrain the ity, and B(λ,T) is the Planck function at temperature T.
distributionofcoremasses,sizesandvelocitydispersionsin Forλ=1.3mm,MAN98adoptκ(1.3mm)=0.005cm2g−1.
Ophiuchus. We use these data to calibrate a multivariate The dust temperature estimates used by MAN98 have
lognormal distribution from which we can draw representa- beenrevisedbyStamatellos,Whitworth&Ward-Thompson
tivecoremasses,sizesandvelocitydispersions.InSection3 (2007). These temperatures, T and T are given
man98 sww07
we describe how we set up the initial conditions for a core; in Table 1. The distance to Ophiuchus used by MAN98,
this includes assigning an ellipsoidal shape to the core and D =160pc,hasalsobeenrevised,byMamajek(2008),
man98
giving it a velocity field and a density profile. In Section to D = 139±6pc. To account for this, each mass given
m08
4 we outline the SPH code used. For each core, we perform by MAN98 must be transformed thus:
simulationswithnoradiativefeedback,continuousradiative
B(1.3mm,T )D2
feedback, and episodic radiative feedback. In Section 5 we M → man98 m08 M , (2.2)
present the results of the simulations and in Section 6 we man98 B(1.3mm,Tsww07)Dm2an98 man98
discuss and summarise the results. Similarly, each size must be transformed thus:
D
R → m08 R . (2.3)
man98 D man98
man98
ThedusttemperaturesT arecoolerthanT ,which
sww07 man98
2 OBSERVATIONS
increases the estimated core masses. However, the revised
Ophiuchus is a star forming region ∼140pc away from the distance Dm08 is less than Dman98, which decreases the core
Sun (Mamajek 2008). It has been well observed in the mil- masses. The combined effect modestly increases the masses
limetre and submillimetre continuum (e.g. Motte, Andre & of cores in Oph-A, Oph-E and Oph-F, whilst leaving the
Neri 1998; Stanke et al. 2006; Simpson, Nutter & Ward-
Thompson 2008), which traces the column density of dust, 1 MAN98 refer to these groups as “cores”, and to the smaller
and allows observers to estimate the masses and sizes of
densercondensationswithinthemas“clumps”.Wehaveadopted
cores.Lada(1999)calculatesthatOphiuchushastotalmass thenownear-universaloppositeconventionthatthesmalldense
Moph∼1to2×103M(cid:12) anddiameterDoph∼1pc.Thecoresin condensationsarecores,andtheyaregroupedinclumps.
c 2013RAS,MNRAS000,1–??
(cid:13)
Simulating star formation in Ophiuchus 3
Region Tman98 Tsww07 Mman98 2.3below,andTable3inonlinematerial)havefreefalltimes
×
(K) (K) t ∼ 4+4 ×104yrs. Thus, whereas interactions between
FF −2
Oph-A 20 11 1.77 coresmayhavebeenimportantindeliveringthedistribution
Oph-B 12 10 1.01 of core properties we observe today, the future evolution of
Oph-C 12 10 1.01
these cores should not be much influenced by interactions.
Oph-D 12 10 1.01
This conclusion is born out in our simulations, where the
Oph-E 15 10 1.40
first protostar usually forms after about one freefall time,
Oph-F 15 10 1.40
andanysubsequentstarformationtendstobeconcentrated
Oph-J - 10 -
inthevicinityofthisprotostar.Asafurtherprecaution,we
Table 1. Dust temperatures from MAN98 and SWW07. The terminate the simulations after 2×105yrs.
fourthcolumngivesthecorrectionfactorappliedtotheMAN98
Treating individual cores in isolation allows us to per-
coremasseswhenweadoptT =Tsww07 andD=139pc. formsimulationswithsufficientlyhighresolutiontocapture
the formation of protostars close to the opacity limit, and
masses of cores in Oph-B, Oph-C and Oph-D essentially to perform many realisations, but it deprives us of the pos-
unchanged (see Table 1). sibility of considering the effects of – for example – accre-
After these adjustments have been made, the MAN98 tionflowsandtidalperturbationsfromthesurroundings.A
CMF peaks at M ∼ 0.3M . This is lower than the large-scalesimulationtreatingthewholeofOphiuchuswith
PEAK (cid:12)
values found for the Ophiuchus cores by other groups, i.e. the same resolution is currently beyond our computational
M ∼0.9M (Stanke et al. 2006) and M ∼0.4M resources,butitshouldbefeasibletosimulateawholeclump
PEAK (cid:12) PEAK (cid:12)
(Simpson,Nutter&Ward-Thompson2008).Thedifferences (e.g. Oph-A).
inM mustbeattributedtoacombinationofthediffer-
PEAK
entwavelengthsused;thedifferentassumptionsmadeabout
dustopacities,dusttemperatures,andsourcedistances;and 2.3 Lognormal distribution
the different procedures used to identify cores, distinguish
Together, MAN98 and ABMP07 provide 61 core masses,
them from the background, and correct for incompleteness
36 mean core radii and 27 core velocity dispersions; 20
at low masses and poor statistics at high masses. We also
coreshaveallthree.Weusethesedatatocalibrateamulti-
notethatstudiesofotherstarformationregions,e.g.Perseus
variate lognormal distribution of core properties. To within
(Enochetal.2008),thePipeNebula(Rathborneetal.2009)
the errors of small number statistics, and modulo the as-
and Aquilla (Ko¨nyves et al. 2010) find CMFs with M
PEAK sumption of a lognormal distribution, this reproduces sta-
between ∼0.6M and ∼1.0M . We have therefore dou-
(cid:12) (cid:12) tistically the observed properties of cores in Ophiuchus,
bled the MAN98 masses, so that the peak of the CMF is
and their correlations. Specifically, we draw values of x ≡
at M ∼0.6M , i.e. between the Stanke et al. (2006)
PEAK (cid:12) (log(M),log(R),log(σnt)) from a probability density
andSimpson,Nutter&Ward-Thompson(2008)values,and
(cid:18) (cid:19)
within the range found for other star formation regions. P(x)= 1 exp −1(x−µ)TΣ−1(x−µ) ,
(2π)3/2|Σ| 2
(2.6)
2.2 ABMP07 where
ABMP07measurethewidthoftheN2H+(1–0)emissionline µM
from 27 prestellar cores in Ophiuchus. Of these 27 cores, µ≡µ  , (2.7)
R
26 have mass estimates from MAN98. The line width gives µ
c
the radial velocity dispersion, σ , of the N H+ molecules
1d 2 and
in the core. In general there are thermal and non-thermal
contributionstothelinewidth,σ2 =σ2+σ2 .Thethermal  σ2 ρ σ σ ρ σ σ 
contribution is given by 1d t nt Σ ≡ ρ σM σ M,Rσ2M R ρM,σnt σMσσnt .
ρ M,RσMσR ρ σR σ R,σnσt2R σnt
σt2 = mkbT , (2.4) M,σnt M σnt R,σnt R σnt σnt (2.8)
mol Hereµ isthearithmeticmeanoflog(x),σ isthestandard
x x
wheremmol =4.8×10−23gisthemassofanN2H+molecule. deviation of log(x) and ρx,y is the Pearson’s correlation co-
The non-thermal contribution, efficient of log(x) and log(y). The values of these terms are
(cid:18) k T (cid:19)1/2 giveninTable2.Fig.1showstheprobabilitydensityfunc-
σ = σ2 − b sww07 , (2.5) tionP(x)projectedthroughallthreeofitsdimensions.Su-
nt 1d m
mol perimposed are the observational data points and 100 data
givesthemagnitudeofmacroscopicgasmotions,i.e.turbu- points drawn randomly from the distribution.
lence, rotation and contraction or expansion. Since we are
adopting the reduced temperatures from SWW07, we de-
rivesomewhathighervaluesofσ thanABMP07,between
nt 3 INITIAL CONDITIONS
∼0.05kms−1 and ∼0.3kms−1.
ABMP07alsomeasurethebulkradialvelocitiesofcores While the mass, size and velocity dispersion of a core are
in Ophiuchus. From the dispersion in the bulk radial veloc- relativelyeasytoquantifyfromobservations,thesamecan-
ities of the cores within a clump, they estimate that, on not be said for the intrinsic three-dimensional shape, den-
average, a core should collide with another core after a few sity profile and internal velocity field. Here we define the
times 105yrs. In contrast, the individual cores (see Section methodology we adopt for assigning these properties.
c 2013RAS,MNRAS000,1–??
(cid:13)
4 O. Lomax, A. P. Whitworth, D. A. Hubber, D. Stamatellos and S. Walch
log(M/M ) log(σnt/kms−1)
⊙
-2.0 -1.5 -1.0 -0.5 0.0 0.5 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2
4.0 4.0
3.8 3.8
3.6 3.6
3.4 3.4
AU) 3.2 3.2 AU)
/ /
R R
g( 3.0 3.0 g(
o o
l l
2.8 2.8
2.6 2.6
2.4 2.4
2.2 2.2
-1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2
-0.4 -0.4 log(σnt/kms−1)
-0.6 -0.6
1)− -0.8 -0.8 1)−
s s
m m
/k -1.0 -1.0 /k
nt nt
σ σ
g( -1.2 -1.2 g(
o o
l l
-1.4 -1.4
-1.6 -1.6
-1.8 -1.8
-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0
log(M/M )
⊙
Figure 1. The multivariate lognormal distribution, P(x) where x=(log(M),log(R),log(σnt)). (a) shows the projections through
log(σnt); (b) through log(M); and (c) through log(R). The concentric ellipses show the 1σ, 2σ and 3σ regions of the distribution.
ThegreencirclesarerandomlydrawnpointsfromP(x).TheredsquaresaretheobservationaldatafromMAN98andABMP07.Because
thereisalimtednumberofcasesforwhichmorethanoneoftheseparametersisknown,thedistributionofobervationalpointsisshifted
relativetothefitteddistribution.
Parameter Value generated with
µM [log(M/M(cid:12))] -0.57
µ [log(R/AU)] 3.11
R
µσnt [log(σnt/kms−1)] -0.95 A=1,
σM [log(M/M(cid:12))] 0.43 B =exp(τGb), (3.1)
σ [log(R/AU)] 0.27
σRσnt [log(σnt/kms−1)] 0.20 C =exp(τGc),
ρ 0.61
M,R
ρM,σnt 0.49 where Gb and Gc are random numbers drawn from a Gaus-
ρ 0.11
R,σnt sian distribution with zero mean and unit standard devia-
Table 2.Arithmeticmeans,standarddeviationsandcorrelation tion.ForOphiuchus,agoodfitisobtainedwithτ=0.6,and
coefficientsoflog(M),log(R)andlog(σnt)forcoresinOphiuchus improvements to the fit obtained by introducing additional
parameters do not appear to be justified (LWC13). Once
the axes have been generated using Eqn. (3.1), they are re-
ordered so that A≥B ≥C and normalised so that A=1.
3.1 Shapes Then, given a mean radius R from the distribution in Eqn.
(2.6), the core axes are set to
Inprojection,thecoresinOphiuchushaveirregular,andof-
tenelongated,shapes,whichshownosignificantcorrelation
with mass, size or velocity dispersion. LWC13 show that, if R
A = ,
one assumes the intrinsic three-dimensional shape of a core core (BC)13
isarandomlyorientedtriaxialellipsoid,theobservedaspect (3.2)
B =BA ,
ratioscanbefittedwithaoneparametermodel.Specifically, core core
C =CA .
therelativesizesoftheprincipalaxesoftheellipsoidcanbe core core
c 2013RAS,MNRAS000,1–??
(cid:13)
Simulating star formation in Ophiuchus 5
3.2 Density profile (k ,k ,k ), in the range 0 to 64. For a power spectrum
x y z
P ∝ k−α the energy density in k-space is E ∝ k−(α+2).
Observationssuggest(e.g.Alves,Lada&Lada2001;Harvey k k
Hence, with α=2, E ∝ k−4 and the energy is strongly
etal.2001;Ladaetal.2008)thatdensecoresoftenapproxi- k
concentratedtowardslowspatialfrequencies(i.e.longwave-
matetothedensityprofileofacriticalBonnor-EbertSphere
lengths).Theamplitude,a,andphase,ϕ,ofeachindividual
(Bonnor 1956), i.e. ρ = ρ e−ψ(ξ), where ρ is the central
C C mode, k, are set by generating a vector G of three random
density,ψ istheIsothermalFunction,ξ isthedimensionless
numbers drawn from a Gaussian distribution having zero
radius, and the boundary is at ξ = 6.451. The FWHM of
B mean and unit variance, and a vector R of three random
the column-density through a critical Bonnor-Ebert Sphere
numbers drawn from a uniform distribution between zero
corresponds to ξ =2.424, so the density at (x,y,z) is
FWHM and one. Then we put
given by
(cid:112)
a(k) = E(k)G, (3.6)
(cid:18) x2 y2 z2 (cid:19)1/2
ξ = ξ + + , (3.3) ϕ(k) = 2πR, (3.7)
FWHM A2 B2 C2
core core core
M ξ e−ψ(ξ) The corresponding contribution to the dimensionless veloc-
ρ = core B , ξ<ξ , (3.4) ity field is given by
4πA B C ψ(cid:48)(ξ ) B
core core core B
vˆ(k) = a(k)cos(ϕ(k))+ia(k)sin(ϕ(k)), (3.8)
where ψ(cid:48) is the first derivative of ψ.
Using this density profile implies that the cores are in v(x) = 1 Re(cid:18)(cid:90) vˆ(k)eikxd3k(cid:19) ; (3.9)
hydrostatic equilibrium, contained by an external pressure. (2π)2
For the simulated cores is this paper, this is generally not
Eqn. (3.9) is evaluated numerically using the fast Fourier
thecase.However,thedensityprofiledoesreproducetheflat
transformlibraryFFTW(Frigo&Johnson2005)toproduce
centralregionandpower-lawenvelopethatisoftenobserved values of the velocity field on a 1283 grid.
in cores (e.g. Kirk, Ward-Thompson & Andr´e 2005; Roy
et al. 2013).
3.3.2 Ordered motions
We introduce large scale rotation and radial motion
3.3 Velocity fields
into the velocity field by simply modifying the phases
Spectroscopic observations of prestellar cores indicate that and symmetries of the k = 1 modes, i.e. k =
theirvelocityfieldsmaycontainorderedrotation(e.g.Good- (1,0,0), (0,1,0)and(0,0,1). Specifically, we require that
man et al. 1993), ordered inward and/or outward motion
   
a(1,0,0) r ω −ω
(e.g. Keto et al. 2006) and random turbulence (e.g. Dubin- x z y
ski, Narayan & Phillips 1995). If the corresponding contri- a(0,1,0)=−ωz ry ωx  ,
a(0,0,1) ω −ω r
butions to the radial velocity dispersion are σ , σ and y x z
rot rad (3.10)
   
σ , the total is given by ϕ(1,0,0) π/2 π/2 π/2
turb
σ2 =σ2 +σ2 +σ2 . (3.5) ϕ(0,1,0)=π/2 π/2 π/2 .
nt rot rad turb ϕ(0,0,1) π/2 π/2 π/2
For the simplest case of a spherically symmetric core ro-
Here (r ,r ,r ) are radial velocity amplitudes and
x y z
tating about an axis perpendicular to the line of sight, it
(ω ,ω ,ω ) are rotational velocity amplitudes; the diame-
x y z
requires two parameters to specify the relative contribu-
ter of the core is scaled to π, so that it corresponds to the
tions to σ from rotation, isotropic pulsation and isotropic
nt central half wavelength of these modified modes. The out-
turbulence. An additional parameter is required if the ro-
come of applying this modification is shown in Fig. 2.
tation axis is not perpendicular to the line of sight. The
The radial velocity amplitudes are proportional to the
problembecomesevenmorecomplicatedifthecoreisellip-
rates at which the core expands along the corresponding
soidal,sincethepulsationsandturbulencearethenunlikely
axis. For example, if r < 1, then the core is contracting
x
tobeisotropic;forexample,contractionislikelytobefastest
along the x-axis, and if r >1, the core is expanding along
x
along the shortest axis (Lin, Mestel & Shu 1965). We have
the x-axis. On the assumption that the core is contracting
no way of estimating specific values for any of these pa-
fastest along its shortest axis, and slowest along its longest
rameters. To overcome this problem, we generate a random
axis, we re-order the amplitudes so that r ≤ r ≤ r , and
x y z
turbulentvelocityfield,characterisedonlybythepowerex-
orient the core so that its ellipsoidal axes, A≥B ≥C, are
ponent (we adopt α=2; cf. Burkert & Bodenheimer 2000)
respectivelyalignedwiththex,yandzaxes.Therotational
and the ratio of compressive to solenoidal modes (we adopt
amplitudesdefinetheangularvelocityω=(ω ,ω ,ω ),giv-
x y z
thethermalratio,i.e.M=1/2),butshiftthecentreofthe
ing the axis and rate of core rotation; the axis of rotation
coresothatitisatthecentreofthelargestmode;theenergy
has a random direction. The values of r , r , r , ω , ω
x y z x y
inthelargestmodeistheninvestedinorderedrotationand
and ω are independently drawn from a Gaussian distribu-
z
compression, whilst the rest is in smaller-scale turbulence.
tion having zero mean and unit variance. This ensures that
the energy in these large-scale components of the velocity
field still subscribes to the same P ∝ k−2 power law as
3.3.1 Turbulent motions k
the smaller-scale ones that deliver internal turbulence. We
To construct a dimensionless turbulent velocity field, we stressthatonlythek=1modesaremodified.Allmodeswith
work in a (2π)3 cubic domain, and populate all modes k>1aregeneratedasdescribedinSection3.3.1.Hencewe
for which the wave-vector, k has integer components, have a velocity field composed of ordered rotation, ordered
c 2013RAS,MNRAS000,1–??
(cid:13)
6 O. Lomax, A. P. Whitworth, D. A. Hubber, D. Stamatellos and S. Walch
(a) (b) (c) (d)
Figure 2. 2D velocity fields generated in an analogous manner to the 3D velocity fields described in Section 3.3.1 and 3.3.2. All
four frames show the region π x π and π y π. From left to right: (a) a compressive velocity field with a(1,0)=( 1,0),
−2 ≤ ≤ 2 −2 ≤ ≤ 2 −
a(0,1)=(0, 1),ϕ(1,0)=(π,0)andϕ(0,1)=(0,π);(b)asolenoidalvelocityfieldwitha(1,0)=(0,1),a(0,1)=( 1,0),ϕ(1,0)=(0,π)
− −
and ϕ(0,1) = (π,0); (c) a turbulent field with modes generated as described in Section 3.3.1, but without the k=1 modes; (d) the
additionofallthreefieldsgivingasuperpositionoflarge-scaleorderedmotionandsmallerturbulentmodes.
radialmotionsand stochasticturbulentmotions,andbased the transport of cooling radiation against dust opacity and
on just two parameters, α and M. molecular opacity, the sublimation of dust, the ionisation
and dissociation of hydrogen and the ionisation of helium.
It has been extensively tested, and reproduces very accu-
3.3.3 Interpolation rately the detailed results of Masunaga & Inutsuka (2000).
To compute the initial velocity of an SPH particle, i, we
assign it a position vector, 4.2 Accretion luminosity
(cid:18) (cid:19)
s = π xi , yi , zi , (3.11) Thedominantcontributiontotheluminosityofaprotostar
i 2 Acore Bcore Ccore is from accretion,
and estimate its dimensionless velocity by interpolation on fGM M˙
L (cid:39) (cid:63) (cid:63) , (4.1)
thegridofcomputedvalues.Finallythedimensionlessveloc- (cid:63) R
(cid:63)
itiesoftheSPHparticlesarescaledsothattheir1Dvelocity
where, f = 0.75 is the fraction of the accreted material’s
dispersion matches a value of σ drawn from Eqn. (2.6).
nt gravitationalenergythatisradiatedfromthesurfaceofthe
protostar(therestispresumedtoberemovedbyradiation,
andbybipolarjetsandoutflows; Offneretal.2009);M is
(cid:63)
4 NUMERICAL METHOD the mass of the protostar; M˙ is the rate of accretion onto
(cid:63)
the protostar; and R = 3R is the approximate radius of
(cid:63) (cid:12)
4.1 SPH code
a protostar (Palla & Stahler 1993).
We use the seren (Hubber et al. 2011) implementation In the sequel, each set of initial conditions is used to
of grad-h SPH (Springel & Hernquist 2002; Price & Mon- run three simulations. In the first, radiative feedback from
aghan 2004) with smoothing parameter η=1.2 (hence typi- sinksisneglected(NRF).Inthesecond,materialassimilated
cally∼56neighbours).GravityiscalculatedwithaBarnes- byasinkispresumedtobeaccretedinstantaneouslybythe
Hut tree and a gadget-style multipole acceptance crite- central protostar, and therefore delivers continuous radia-
rion (Springel & Hernquist 2002). The code invokes mul- tive feedback (CRF). In the third, material assimilated by a
tiple particle timesteps and a leapfrog drift-kick-drift inte- sinkisassumedtoaccumulateinaninneraccretiondisc,and
grationscheme.Time-dependentartificialviscosityisimple- thenaccreteontothecentralstarinaburst,givingepisodic
mented using the Morris & Monaghan (1997) formulation radiativefeedback(ERF).TheNRFsimulationsareperformed
with α=1, α =0.1 and β =2α. purelyforcomparison,andourmainconcerniswiththedif-
min
All simulations are set up with m = 10−5M , so ferencebetweencontinuousandepisodicradiativefeedback.
sph (cid:12)
the mass resolution is ∼ 10−3M . A gas concentration In order for the disc around a protostar to fragment
(cid:12)
that is gravitationally bound and has density higher than gravitationally, at radius R, two criteria must be fulfilled.
ρ = 10−9gcm−3 is replaced with a sink particle (Bate, First, gravity must be able to overcome thermal and cen-
sink
Bonnell & Price 1995) having radius R ≈ 0.2AU. Any trifugal support,
sink
SPH particles that subsequently enter this radius and are c(R)κ(R)
gravitationally bound to the sink are assimilated by it. We Σ(R) >∼ πG , (4.2)
have not used the more sophisticated NewSinks method,
where Σ(R) is the local surface-density, c(R) the sound
presentedbyHubber,Walch&Whitworth(2013),however,
speed,andκ(R)theepicyclicfrequency;(Toomre1964).Sec-
for this value of ρ , the two methods are reasonably well
sink ond,thecondensingfragmentmustcoolfastenoughtoavoid
converged.
undergoing an adiabatic bounce and then being sheared
Heating and cooling of the gas are treated with the
apart, i.e.
schemedescribedinStamatellosetal.(2007).Thistakesac-
count of heating by cosmic rays and background radiation, t < t , (4.3)
cool ∼ orb
c 2013RAS,MNRAS000,1–??
(cid:13)
Simulating star formation in Ophiuchus 7
where t is the gas cooling time, and t is the orbital
cool orb
period (Gammie 2001). Simulations that include instanta-
neousaccretionandcontinuousradiativefeedback(e.g.Bate 1 NRF (a)
2009b; Krumholz 2006; Krumholz et al. 2010; Offner et al.
2009, 2010; Urban, Martel & Evans 2010) find that proto- ))
M
sthteelrleafrordeisdcos naoret sspoawwnarlmowt-hmaatssthHe-ybucranninnogtsftraargsmoernbtroawndn (log( 0.1
P
dwarfs.
However, there is observational evidence suggesting
that accretion onto a protostar may sometimes be episodic, 0.01
1 ERF (b)
i.e. it may occur intermittently, in short, intense bursts
(e.g. Herbig 1977; Dopita 1978; Reipurth 1989; Hartmann
& Kenyon 1996; Greene, Aspin & Reipurth 2008). In par- ))
M
ticular, FU Ori-type stars (e.g. Herbig 1977; Greene, Aspin og( 0.1
& Reipurth 2008; Peneva et al. 2010; Green et al. 2011; P(l
Principe et al. 2013) exhibit large increases in their bright-
ness that last for between a few and a few hundred years.
The physical cause of episodic accretion is uncertain. The 0.01
1 CRF (c)
critical issue is the duty-cycle, specifically the duration of
the low-luminosity phase: if this is longer than the cooling
timescale, typically a few times 103yrs, the disc can cool M))
and fragment. og( 0.1
(l
For the purpose of the simulations presented here P
we adopt the phenomenological ERF model of Stamatel-
los, Whitworth & Hubber (2011, 2012) (hereafter SWH11),
0.01
which is based on the calculations by Zhu, Hartmann & 1 NRF (d)
ERF
Gammie (2009, 2010) and Zhu et al. (2010). In the outer CRF
disc(outisidethesink),angularmomentumisredistributed
))
bygravitationaltorques,allowingmaterialtospiralinwards M
andintothesink.However,insidethesinkthematerialnear (log( 0.1
P
the protostar is so hot that the inner disc becomes gravita-
tionally stable, and there are no gravitational torques to
transport angular momentum outwards. Consequently this
0.01
material cannot spiral inwards any further, and instead ac- -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0
log(M/M )
cumulatesintheinnerdisc.Thiscontinuesuntiltheaccumu- (cid:12)
lated material becomes so hot that it is ionised thermally,
sufficiently to couple to the magnetic field. The magneto-
rotational instability (MRI) then cuts in and the resulting Figure3.Protostarmassfunctions.Thetopthreeframes(a,b,c)
torquesenablematerialtospiralinwardsandontothestar, givethemassfunctionsfromsimulationswith,respectively,NRF,
giving a short burst of very high luminosity. The downtime ERF and CRF. The black histograms have bins equally spaced in
– during which material is accumulating in the inner disc, log(M) and the red lines are kernel smoothed density functions.
the luminosity is low, and the outer disc may be able to The bottom frame (d) shows the smoothed NRF, ERF and CRF
densityfunctionsonasingleplot.
fragment – is given by
(cid:18) M (cid:19)2/3(cid:18) M˙ (cid:19)−8/9
t ∼13kyr (cid:63) sink ; (4.4)
down 0.2M 10−5M yr−1
(cid:12) (cid:12)
M˙ is the rate at which material flows into the sink.
sink
5.1 Protostellar masses
5 RESULTS
Fig.3shows,ashistogramsandsmoothedprobabilitydensi-
Eachcoreisevolvedfor2×105yrs,whichislongerthanthe ties2, the protostellar mass functions from the simulations.
freefall time of the least-dense core. Of the 100 cores sam- Both the NRF and ERF simulations deliver a wide range of
pled from the lognormal distribution, 60 collapsed to form masses from a few Jupiter masses to between one and two
at least one star; note that the simulations with different solar masses. The CRF simulations deliver roughly the same
feedback are identical up until this point. The star forma- uppermasslimit,butnoprotostarsbelowabout50Jupiter
tion efficiencies and numbers of protostars formed per core masses. In the NRF simulations, the mean efficiency (total
aregiveninTable3intheonlinematerial.Inthefigureswe mass of stars formed divided by core mass) is 71%, with on
use the naming convention NUM FBK, where NUM is the core average 3.3 stars and 3.4 brown dwarfs per core. In the ERF
number (identifying the initial conditions) and FBK is the simulations,theefficiencyisagain71%,withonaverage3.1
feedback type. starsand1.5browndwarfspercore.IntheCRFsimulations,
c 2013RAS,MNRAS000,1–??
(cid:13)
8 O. Lomax, A. P. Whitworth, D. A. Hubber, D. Stamatellos and S. Walch
theefficiencyis59%,withonaverage1.6starsand0.1brown
dwarfs per core.
1 NRF (a)
Kroupa(2001)
Chabrier(2005)
Ratio of stars to brown dwarfs ))
M
Andersenetal.(2008)haveperformedasurveyofstarform- (log( 0.1
P
ing regions and find that low mass stars outnumber brown
dwarfs by a factor
N(0.08M <M ≤1.0M ) 0.01
A= (cid:12) (cid:12) =4.3±1.6. (5.3) 1 ERF (b)
N(0.03M <M ≤0.08M ) Kroupa(2001)
(cid:12) (cid:12) Chabrier(2005)
In the simulations presented here, A =2.2±0.3, A =
nrf erf ))
3.9±0.6andA =17±8.ThisconfirmsthatNRFproduces M
toomanybrowcnrfdwarfs,andsuggeststhatCRFproducestoo (log( 0.1
P
few. ERF appears to produce about the right number.
0.01
Comparison to the IMF 1 CRF (c)
Kroupa(2001)
Chabrier(2005)
InFig.4,wehaveplottedthemassdistributionsfromFig.3
with the Kroupa (2001; K01) and Chabrier (2005; C05) fits ))
M
ttohethCe0I5MIFM.FTahteMK0pCe10a5kIM=F0p.2eaMk(cid:12)s.atInMtphKe0ea1ks=imu0.l0a8tiMon(cid:12)s,wainthd P(log( 0.1
NRF, ERF and CRF, the mass distributions peak at Mnrf ≈
peak
0.1M , Merf ≈0.2M and Mcrf ≈0.4M , respectively.
(cid:12) peak (cid:12) peak (cid:12)
ThepeaksoftheNRFandERFmassdistributionsareingood 0.01
-3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0
agreementwiththoseofK01andC05respectively;thepeak
log(M/M )
of the CRF mass distribution appears to be rather too high. (cid:12)
There are no high-mass cores in Ophiuchus, and therefore
none of the simulated IMFs match the observed IMFs at Figure 4. As Fig. 3, but here the simulated protostellar mass
high masses. functions (in black) are compared with the K01 and C05 IMFs
(inredandbluerespectively).Thethreeframesshowprotostellar
massfunctionforsimulationswith(a)NRF,(b)ERF,and(c)CRF.
ThedashedpartsoftheK01andC05IMFsareinthestellarmass
Summary: masses
regime (M 0.08M(cid:12)) and the dotted parts are in the brown
≥
For masses less than ∼1M , the ERF simulations deliver dwarf mass regime (M < 0.08M(cid:12)), which is not as well con-
(cid:12)
strainedasthestellarregime.Allmassfunctionsarenormalised
a protostellar mass function which is compatible with ob- sothat(cid:82)+∞P(log(M))dlog(M)=1.
servedIMFs.Incontrast,theCRFsimulationsappeartopro- −∞
ducetoofewlow-massH-burningstarsandbrowndwarfs.In
the ERF simulations each core produces on average N ≈
ERF 5.2 Fragmentation
4.6 stars with a star formation efficiency of η ≈ 0.7.
ERF
ThisisinagreementwiththestatisticalinferenceofHolman Mechanisms
et al. (2013) that values of N =4.3±0.4 and η=1.0±0.3
In many simulations, the core collapses to form a central
are required to reproduce the observed multiplicity statis-
protostarsurroundedbyanaccretiondisc.Ifthediscissuf-
tics (η >1 simply implies that the cores accrete additional
ficiently massive and cold, it fragments to form additional
massastheyformstars).IntheCRFsimulationsthenumber
protostars, typically with masses smaller than the central
of stars formed is significantly lower, N = 1.7, and the
CRF protostar.AnexampleofthisisshowninFig.5.Weseethat
efficiency is a little lower, η =0.59.
CRF thetimescaleofdiscfragmentationisoforder103yrs.This
mechanismiswelldocumented(e.g.Stamatellos,Hubber&
Whitworth 2007; Stamatellos & Whitworth 2008, 2009a,b;
2 Thesmoothedprobabilitydensityisgivenby Thies et al. 2010; Stamatellos et al. 2011) and has been
1 (cid:88)N 1 (cid:18) (log(M) log(Mi))2(cid:19) shown to reproduce the observed properties of stars and
P(log(M))= exp − . brown dwarfs rather well (e.g. Stamatellos & Whitworth
N √2πh2 − 2h2
i=1 2009a).
(5.1)
In most simulations the central regions are fed by fila-
IntheexpectationthatP(log(M))isroughlylognormal,weset
mentary streams of material, and in some cases these fila-
(cid:18)4σˆ5(cid:19)51 ments fragment before a central protostar forms; this also
h= (5.2)
3N occurs on time scales of order 103yrs. An example of this
is shown in Fig. 6. This core has a relatively high mass
where σˆ is the standard deviation of log(M) (Silverman 1998).
Thisisausefulmethodofextractingfromnoisyhistogramslarge (M=3.3M(cid:12)) and a low ratio of turbulent to gravitational
scalefeaturesthatareotherwisehardtosee. energy(α =0.07).Thecorecollapsesintoafilamentbe-
turb
c 2013RAS,MNRAS000,1–??
(cid:13)
Simulating star formation in Ophiuchus 9
Figure 5. Column density plots from simulation 002 ERF show- Figure 7. Column density plots from the 022 XXX simulations.
ing disc fragmentation. From left to right, then top to bottom, The left-hand column gives snapshots from 022 NRF, the middle
the times are t=17,18,19and20kyr. This core fragments into column from 022 ERF, and the right-hand column from 022 CRF.
three stars and three brown dwarfs. This figure and subsequent Fromtoptobottom,thetimesaret=15,18and23kyr.022 NRF
simulationsnapshotsarerenderedusingSplash(Price2007). formsfourstarsandninebrowndwarfs,022 ERFformsfourstars,
and022 CRFformsasinglestar.
simulations. Second, in the NRF simulations there is a sec-
ondpeakaround∼0.01M .Themainpeaksarepopulated
(cid:12)
by protostars that form from core collapse. The secondary
peak in the NRF simulations arises from prolific disc frag-
mentation.Thisproducesanabundanceoflowmassobjects,
manyofwhichareejectedfromthecore(e.g.Bate,Bonnell
& Bromm 2002). Once ejected, accretion ceases and their
masses are frozen at low values.
In the ERF simulations, high-luminosity outbursts peri-
odicallyraisethetemperatureofthegassurroundingapro-
tostar,stabilisingitagainstgravitationalfragmentation.Ini-
tially,theseoutburstsoccuratintervalsof∆t ∼104yrs,
down
whereasgravitationalinstabilitiesdevelopontimesscalesof
∆t ∼ 103yrs; therefore there are windows of opportunity
gi
for fragmentation between outbursts. However, as further
protostars form, the intervals between outbursts decrease
until eventually ∆t <∆t . Thereafter, further star for-
down∼ gi
mation is inhibited by regular luminosity outbursts. Thus
Figure6.Columndensityplotsfromsimulation098 ERFshowing
ERF allows fragmentation to occur, but then regulates the
the filament fragmentation. From top to bottom, the times are
process after a few protostars have formed. Because fewer
t = 19,23,27and36kyrs. This core fragments into eight stars
low-mass protostars are formed, there are fewer dynamical
andtwobrowndwarfs.
ejections,andmoremassisavailabletoaccreteontothecen-
tral protostar (cf. Girichidis et al. 2012), so the main peak
causeitisparticularlyelongated,withitslongestaxismore is enhanced and the secondary peak is removed.
than twice the length of the other two. The filament then In the CRF simulations, once the first protostar has
spawnstenprotostars.ThisconfirmstheresultsofGoodwin, formed, the remaining gas is continuously heated and very
Whitworth&Ward-Thompson(2004),whoshowedthatlow littlefurtherfragmentationoccurs,sothesesimulationsusu-
levelsofturbulencearesufficienttofragmentsphericalcores. allyformsingleprotostarsorbinaries.Fig.7showsanexam-
pleofacoreevolvingwithNRF,ERFandCRF;theusualtrend
isthatthemostprotostarsareformedwithNRF,followedby
Feedback effects ERF, then CRF.
Fig.7showsacomparisonofcolumndensityplotsfromsim-
ulationsusingdifferentfeedbackmechanisms.Therearetwo
Number of objects formed
maindifferencesbetweentheprotostellarmassdistributions
fromtheNRFandERFsimulations.First,themainpeakshifts Fig. 8 shows that the number of protostars formed per
from∼0.1M intheNRFsimulations,to∼0.2M intheERF core is correlated with the core mass. Cores with masses
(cid:12) (cid:12)
c 2013RAS,MNRAS000,1–??
(cid:13)
10 O. Lomax, A. P. Whitworth, D. A. Hubber, D. Stamatellos and S. Walch
0.0
NRF NRF
10.0 (a) -0.5
E )(cid:12)
OR M
/NC M/? -1.0
N? 1.0 og(
l
-1.5
(a)
0.1 -2.0
ERF ERF
10.0 (b) -0.5
E )(cid:12)
OR M
/NC M/? -1.0
N? 1.0 og(
l
-1.5
(b)
0.1 -2.0
CRF CRF
10.0 (c) -0.5
E )(cid:12)
OR M
/NC M/? -1.0
N? 1.0 og(
l
-1.5
(c)
0.1 -2.0
-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0
log(MCORE/M ) log(MCORE/M )
(cid:12) (cid:12)
Figure 8. The number of protostars formed per core as a func- Figure 9. The distribution of protostar masses as a function
tion of core mass. Error bars reflect the Poisson uncertainties in of core mass. The solid line gives the geometric mean mass and
the numbers of cores and sinks, per bin. From top to bottom, standard error. The shaded area shows the geometric standard
simulationsinvoking(a)NRF,(b)ERFand(c)CRF. deviation. The dashed line gives the arithmetic mean protostar
mass.Fromtoptobottom,simulationsinvoking(a)NRF,(b)ERF
and(c)CRF.
Fragmentation: summary
between 0.1M and 0.3M typically form single proto-
(cid:12) (cid:12)
AsnotedbyStamatellos,Whitworth&Hubber(2012),ERF
stars, whereas cores with masses between 3M and 10M
(cid:12) (cid:12)
providesawindowofopportunityfordiscfragmentation.As
typically produce more than ten protostars with ERF, and
further protostars form, they too radiate and the window
about five with CRF. For convenience, we will refer to cores
eventually becomes too short to allow fragmentation. This
with M ≤ 1M as low-mass cores, and cores with
core (cid:12) meansthat,unlikeCRF,ERFpermitsbrowndwarfstoformby
M >1M as high-mass cores.
core (cid:12) disc fragmentation, but, unlike NRF, brown dwarf formation
Fig. 9 shows the distribution of protostar masses as a
is regulated to a level that agrees better with observation.
function of core mass. In the NRF and ERF simulations, the
Irrespectiveofthetypeoffeedbackinvoked,thesesimu-
arithmetic mean mass is roughly constant over the entire
range of core masses, viz. M¯nrf = 0.12M and M¯erf = lationsdonotsupportthenotionthattheshapeoftheIMF
(cid:63) (cid:12) (cid:63) is inherited from the shape of the CMF. The mean mass
0.19M .AlsoshownonFig.9isthegeometricmeanofthe
(cid:12) of a protostar formed in a low-mass core is approximately
protostarmassesanditsstandarddeviation.Thisshowsthat
the same as one formed in a high-mass core. The difference
high-masscoresformawiderrangeofmassesthanlow-mass
is that high-mass cores tend to produce more protostars,
cores. In low-mass cores, there is little spread in the proto-
roughly in proportion to their mass.
star masses because most of the available gas goes into the
Competitive accretion occurs in our simulations, but
central protostar.A disc may form,but itis oftennot mas-
only on the scale of individual cores having sufficient mass
siveenoughtofragment.High-masscoreshavesufficientgas
(M >M ) to form multiple protostars.
to form more massive protostars, and there is often enough CORE (cid:12)
gas remaining to build a massive disc, which can then frag-
mentintolow-massprotostars.Finalmassesaredetermined
5.3 Limitations
by competitive accretion for the remaing core mass.
In simulations with CRF, low-mass cores tend to form SincethemostmassivecoreinOphiuchushasmass∼5M
(cid:12)
singleprotostars,andhigh-masscorestendtoformmultiple (and the next most massive ∼ 3M ), and since we evolve
(cid:12)
protostars, but in smaller numbers than with ERF. our cores in isolation (i.e. they do not accrete extra mate-
c 2013RAS,MNRAS000,1–??
(cid:13)