Star Cluster Formation from Turbulent Clumps. I. The Fast Formation Limit

8.2 MB·

by  Juan P. Farias

Draft version March 3, 2017
PreprinttypesetusingLATEXstyleemulateapjv.01/23/15
STAR CLUSTER FORMATION FROM TURBULENT CLUMPS.
I. THE FAST FORMATION LIMIT
Juan P. Farias1,(cid:63), Jonathan C. Tan1 and Sourav Chatterjee2
1AstronomyDepartment,UniversityofFlorida,Gainesville,FL32611,USA
2CenterforInterdisciplinaryExploration&ResearchinAstrophysics(CIERA)
Physics&Astronomy,NorthwesternUniversity,Evanston,IL60202,USA
Draft version March 3, 2017
ABSTRACT
7
We investigate the formation and early evolution of star clusters assuming that they form from a
1
turbulent starless clump of given mass bounded inside a parent self-gravitating molecular cloud char-
0
acterized by a particular mass surface density. As a first step we assume instantaneous star cluster
2
formation and gas expulsion. We draw our initial conditions from observed properties of starless

r clumps. We follow the early evolution of the clusters up to 20 Myr, investigating effects of different
a
starformationefficiencies,primordialbinaryfractionsandeccentricitiesandprimordialmasssegrega-
M
tion levels. We investigate clumps with initial masses of M =3000M embedded in ambient cloud
cl (cid:12)
environments with mass surface densities, Σ = 0.1 and 1gcm−2. We show that these models of
1 cloud
faststarclusterformationresult,inthefiducialcase,inclustersthatexpandrapidly,evenconsidering

only the bound members. Clusters formed from higher Σ environments tend to expand more
] cloud
A quickly, so are soon larger than clusters born from lower Σ conditions. To form a young cluster
cloud
ofagivenage, stellarmassandmasssurfacedensity, thesemodelsneedtoassumeaparentmolecular
G
clump that is many times denser, which is unrealistic compared to observed systems. We also show
h. that in these models the initial binary properties are only slightly modified by interactions, meaning
p that binary properties, e.g., at 20 Myr, are very similar to those at birth. With this study we set up
- the basis of future work where we will investigate more realistic models of star formation compared
o to this instantaneous, baseline case.
r
t Keywords: methods: numerical — galaxies: star formation, star clusters
s
a
[
1. INTRODUCTION from more massive stars (e.g., Dale et al. 2015); contin-

2 Most stars tend to form together in clusters (e.g., ued infall of gas to the clump; dispersal of clump gas by
v feedback;dynamicalinteractionsamongtheformingand
Gutermuth et al. 2009), which are created from over
1 recently formed stars as they orbit in the protocluster
dense gas clumps, typically found in giant molecular
0 potential (e.g., Chatterjee & Tan 2012). This is a com-
clouds(GMCs)(e.g.,McKee&Ostriker2007). Thusun-
7 plicated, multi-scale problem, the full solution of which
derstandinghowstarsformcomeswiththedirectneedto
0 is beyond current computational capabilities. Thus ap-
understand how and where star clusters form. It is still
0 proximate models are necessary. By investigating how
a matter of debate if star cluster formation depends on
1. properties of the GMC environment or not. While theo- the outcome of star cluster formation depends on the
0 retical studies have taught us the essential physical pro- adopted approximations we can learn which processes
7 cesses that determine a star cluster’s evolution after gas are most important.
1 Our approach in this paper and subsequent papers in
is dispersed, the transition from the dense star-forming
: this series is to follow the dynamics of formed stars,
v clump to the star cluster that emerges from the gas is
including binary properties, accurately via direct N-
i not yet well understood (see, e.g., Banerjee & Kroupa
X body integrations, but approximate various models for
2015, for a review). In particular, it is debated whether
how individual stars are born within the star-forming
r the process is slow with the clump evolving in a quasi-
a clump. OurinitialconditionsarebasedontheTurbulent
equilibriumstate(Tanetal.2006;Nakamura&Li2007),
Core/Clump model of McKee & Tan (2003) (hereafter
orveryrapidwithstarclusterformationoccurringinjust
MT03), which approximates clumps as singular poly-
a crossing time of the system (Elmegreen 2000, 2007;
tropic virialized spheres that are in pressure equilibrium
Hartmann & Burkert 2007).
with their surrounding cloud medium. This surrounding
There are numerous physical processes potentially in-
cloudisalsoassumedtobeself-gravitatingsoitsambient
volved in such a transition from gas clump to star clus-
pressure is P¯ ∼GΣ2, where Σ is the cloud mass surface
ter, including: fragmentation of a magnetized, turbulent
density—the main variable describing different environ-
andself-gravitatingmediumtoapopulationofpre-stellar
mental conditions.
cores; collapse of these cores via rotationally supported
In this first paper, we start with the simplest approx-
disks to single or multiple star systems; feedback from
imation for star formation, i.e., instantaneous formation
the forming stars, especially protostellar outflows that
ofthestellarpopulationfromtheinitialgasclumpalong
can maintain turbulence in the clump (Nakamura & Li
with simultaneous, instantaneous expulsion of the re-
2007, 2014), and eventually radiative feedback processes
maining gas that is not incorporated into stars. This
(cid:63)E-mail:jfarias@ufl.edu approximation has often been adopted by previous N-
2 Farias, Tan & Chatterjee
bodystudies(e.g.,Bastian&Goodwin2006;Parkeretal. 2. PRE-CLUSTERCLUMPSASINITIALCONDITIONS
2014; Pfalzner et al. 2015). However, in comparison
The initial conditions for star clusters are constrained
to these previous studies our work is distinguished by
by the observed properties of dense gas clumps within
(1) adopting initial conditions that have been explicitly
GMCs. These locations are also expected to be the
developed for self-gravitating gas clumps (i.e., singular
sites of future massive star formation. The properties of
polytropicspheresasapproximationsforturbulent,mag-
theseclumpshavebeensummarizedbyTanetal.(2014),
netizedclumps);(2)followingthefullevolutionofbinary
based on Galactic observational studies of Infrared Dark
systems.
Clouds (IRDCs) (e.g., Rathborne et al. 2006; Butler &
A number of authors have studied the dynamics of
Tan2009,2012)andmm/sub-mmdustcontinuumemis-
binaries in star clusters using numerical models (e.g.,
sion and molecular line surveys of clumps (e.g., Schuller
Kroupaetal.1999;Kroupa&Burkert2001;Parkeretal.
et al. 2009; Ginsburg et al. 2012; Ma et al. 2013). There
2009;Kouwenhovenetal.2010;Parkeretal.2011;Kacz-
are a range of clump masses observed from ∼100M to
marek et al. 2011), focusing on various aspects of their ∼105M . (cid:12)
dynamics. In our work we follow the evolution of binary (cid:12)
In the fiducial Turbulent Clump model, the mass sur-
properties due to stellar interactions and stellar evolu-
face density of the clump of interest is only a factor
tion in clusters forming in different environments. We
of 1.22 times higher than that of its surrounding cloud
use and will use these results to constrain assumptions
(MT03). Therefore the observed surface densities of
aboutthestarclusterprocess,e.g.,itsduration,andindi-
clumps give a reasonable estimate of the mass surface
vidual star formation processes, e.g., how binaries form,
densitiesoftheambientclouds,Σ ,whichistheother
i.e., the primordial binary properties. We also examine cloud
main parameter needed to set up the initial conditions
the role of binaries in producing dynamical and binary
of the models. The observed values of mass surface den-
supernova ejections, especially fast runaway stars.
sityofGalacticclumpsandprotoclustersareintherange
We anticipate that, in reality, star cluster formation from ∼0.03–1gcm−2.
is likely to proceed in a relatively gradual manner, i.e.,
In the setup of our initial conditions in this paper we
taking at least several and perhaps many local free-fall
makeseveralsimplifyingassumptionsasfirststepsinde-
timescales of the gas clump. Modeling this process of
scribing the complexity of star cluster formation: i) the
gradual star cluster formation will be the topic of the
parentclumpisisolated(i.e., noexternaltidalfields); ii)
second paper in this series. One aspect that is needed in
the clump is in hydrostatic and virial equilibrium with
such a model is the gradual evolution of the potential of
thestructureofasingularpolytropicsphere(MT03);iii)
the natal gas clump, which can be approximated via a
stars are born with the same velocities as their parent
simpleanalyticrelation. Suchanapproachofanevolving
gas,sothattheirvelocitydispersionprofileisthesameas
approximatebackgroundgaspotentialhasbeenadopted
that of the initial gas; iv) all stars form instantaneously
previouslyby(e.g.,Tutukov1978;Ladaetal.1984;Geyer
andtheremaininggasisalsoexpelledinstantaneouslyat
&Burkert2001;Boily&Kroupa2003;Chen&Ko2009;
this time; iv) the star formation efficiency (SFE) is spa-
Smith et al. 2011; Farias et al. 2015; Brinkmann et al.
tially constant, which means that stars follow the same
2016). However, in all these studies the full stellar pop-
spatialdistributionastheinitialgas;v)thereisnoinitial
ulation was introduced instantaneously at the beginning
spatial or kinematic substructure given to the stars, ex-
of the simulation. None of these studies included a full
cept that which results from random, Poisson sampling;
treatmentofbinaries,i.e.,includingasignificantfraction
vi) following an initial test model of equal mass stars,
of primordial binaries.
a standard Kroupa IMF for the stars (Kroupa 2001)
Thesecondaspectisthegradualformationofthestel-
is adopted, then with various binary properties inves-
larpopulationinthecluster. Thiswillbethemainfocus
tigated. Itshouldberememberedthatthesearestarting
of Paper II in our series. There have so far been rela-
assumptions and that many of these will be relaxed in
tively few studies using such an approach. Proszkow &
subsequentinvestigations. However,firstthebehaviorof
Adams (2009) presented a study involving gradual star
this simplest, idealized model needs to be understood.
formation and then the early evolution of young clus-
ters, focusing on low to intermediate mass clusters (100
2.1. Initial stellar phase-space distributions
to 3000 members), without primordial binaries. How-
ever, the stellar densities of their models were relatively First we define the physical and kinematic properties
low, limitingtheeffectsofstellarinteractions. Theyalso of the pre-cluster clump, i.e., mass, size, density pro-
did not include stellar evolution and only simulated up file and velocity dispersion profile. Stars are born from
to 10 Myr. thisclumpandinitiallyfollowthesamephase-space(po-
These methods of N-body modeling can be con- sition, velocity) distribution. MT03 characterize pre-
trasted with other approaches to simulating star cluster clusterclumpsandpre-stellarcoresassingularpolytropic
formation, e.g., simulations that follow the (magneto- spheres in virial and hydrostatic equilibrium. The den-
)hydrodynamics of the collapsing clump (e.g., Price & sity profile of such clumps is then:
Bate2009;Padoanetal.2012;Myersetal.2014;Padoan
et al. 2014). Such simulations must still implement sub- (cid:18) r (cid:19)−kρ
ρ (r)=ρ , (1)
grid models for how stars form and inject feedback into cl s,cl R
cl
the gas. Typically they do not have the resolution to
accurately follow binary orbital evolution. Still, these whereρ isthedensityatthesurfaceoftheclump,i.e.,
s,cl
are complementary approaches to those based on pure at radius R . We adopt k = 1.5 as a fiducial value,
cl ρ
N-body approaches and comparison of the results of the i.e., the same as that of MT03 who made their choice
different methods will be instructive. based on observations of clumps reported at the time.
Star Cluster Formation from Turbulent Clumps 3
No significant difference has been found in later mea- the desired initial density profile, no matter the nature
surements performed by Butler & Tan (2012) in IRDCs of the stellar IMF.
where they found k (cid:39) 1.6. We thus keep k = 1.5 as Thekinematicpropertiesoftheclumparespecifiedby
ρ ρ
ourfiducialvalueforsimplicityandconsistencywiththe theconditionofvirialequilibrium. Inavirializedclump,
previousanalysisofMT03. Thedensityatthesurfaceof the velocity dispersion σ scales in the same manner as
the clump can be expressed as the effective sound speed c ≡ (P/ρ)1/2. Therefore, the
velocity dispersion profile of the clump is:
(3−k )M
ρ = ρ cl, (2)
s,cl 4πR3 (cid:18) r (cid:19)(2−kρ)/2
cl σ (r)=σ , (7)
cl cl,s R
where M =M(r <R ) is the total mass of the clump. cl
cl cl
The radius of a clump in virial equilibrium and pres- where σ is the velocity dispersion at the surface of
cl,s
sure equilibrium with its surroundings, i.e., a larger self- the clump. The presence of large-scale magnetic fields
gravitatingcloudofgivenmasssurfacedensity,Σ ,is canprovidesomesupporttotheclumpsothatasmaller
cloud
(MT03; Tan et al. 2013, hereafter T13) turbulent velocity dispersion is needed to achieve virial
equilibrium. Theeffectofmagneticfieldsonthestability
(cid:18) A (cid:19)1/4(cid:18) M (cid:19)1/2 of the clump can be expressed as φ ≡(cid:104)c2(cid:105)/(cid:104)σ2(cid:105). Then,
R =0.50 cl B
cl kpφP,clφP¯ 3000M(cid:12) the velocity dispersion at the surface is:
×(cid:18)1Σgcclomud−2(cid:19)−1/2 pc (3) σcl,s=5.08(cid:18)φAPk,P2clφφ4BP¯(cid:19)1/8(cid:18)300M0cMl (cid:12)(cid:19)1/4
→0.37Mc1l/,32000Σ−clo1u/d2,1pc (4) ×(cid:18) Σcloud (cid:19)1/4 km s−1 (8)
1 g cm−2
where k = 2(k −1) is the power law exponent of the
p ρ
pressure (P) within the clump; φP,cl is the ratio be- where we use φB = 2.8 as a fiducial value, which is the
tween the pressure at the surface of the clump (P ) value for regions with an Alfv´en Mach number of 1 (see
s,cl
annodrmtahleizmateioanncpornesstsaunret,i∼nsiOde(1t)h,eincltohuedr,ePl¯actlioound;P¯φP¯ is≡a AprpopfielendoifxtAhe2 sotfaMrsTt0h3e)n. fTolhloewisniEtiqa.l v7e.loTcihtye dinisdpiveirdsiuoanl
cloud stellarvelocitiesarethenassignedvelocitiesineachofthe
φP¯GΣ2cloud; and A=(3−kρ)(kρ−1)fg →3/4, since we x, y and z directions independently by drawing from a
assume the clump is initially starless so f =1. As fidu-
g Gaussiancenteredatzerowithwidthσ(r). Notethatthe
cialvalueswechooseφP,cl =2andφP¯ =1.32(seeMT03; massaveragedvelocitydispersionoftheclump/clusteris
T13). Thus the structural properties of the fiducial
(T13)
clump are specified by two parameters: M and Σ .
cl cloud
aInndthinisvsetsutidgyatweeΣwclioluldfo=cu0s.1onantdhe1cgacsmeo−f2,Mwchli=ch3a0r0e0rMep(cid:12)- σcl=28(3−−3kkρ)σcl,s → 76σcl,s. (9)
resentativeoftherangeofvaluesobservedinrealclumps ρ
(Tan et al. 2014). The resulting velocity distribution has the form of a
√
Thetotalmassofstars,M∗,thatformfromtheclump Maxwell-Boltzmann distribution with σ3D = 3σcl but
is given with a one dimensional velocity dispersion profile as in
equation 7. The properties of the clumps, i.e., the low
M =(cid:15)M , (5)
∗ cl and high Σ cases, are summarized in Table 1.
where (cid:15) is the overall SFE. We will consider a range of
2.2. The primordial binary population
valuesfrom(cid:15)=0.1to1,withafiducialvalueof0.5. The
stars are assumed to have the same structural profile as Observational evidence shows that about half of star
theclump,i.e.,(cid:15)isindependentofradiallocation. Then, systems in the field are binaries (e.g., Duquennoy &
the density profile of the stars is simply: ρ (r)=(cid:15)ρ (r). Mayor 1991; Fischer & Marcy 1992; Mason et al. 1998;
∗ cl
The cumulative radial mass distribution of the stars is Preibisch et al. 1999; Close et al. 2003; Basri & Reiners
2006; Raghavan et al. 2010). Given the densities of star-
(cid:18) r (cid:19)3−kρ forming clumps and young star clusters, it is likely that
M (r <R )=M . (6)
∗ cl ∗ R mostofthesebinarieswhereborntogetherinsideindivid-
cl
ual cores, rather then forming via subsequent dynamical
To set up the positions of the stars, we first create the interactions (e.g., Parker & Meyer 2014). However, this
mass sample from a given IMF (including binary com- isaquestionthatoursimulationswillbeabletoaddress.
panions), i.e., in the fiducial case a standard Kroupa Theoretically,afullunderstandingofbinaryformation
IMF (Kroupa 2001) with individual masses in the range from a collapsing gas core is likely to require a full non-
0.01M <m <100M ,inrandomorder. Next,wecre- ideal MHD treatment to resolve formation of the accre-
(cid:12) i (cid:12)
ateacumulativemassarrayfromtheprevioussampleto tiondiskandthenlateritspotentialfragmentation. The
then choose the individual distance from the center r difficulty of this problem means that essentially the sta-
i
according to Eq. 6 (in the case of 100% mass segrega- tistical properties of primordial binaries are very uncer-
tion,seebelow,massesaresortedfromthemostmassive tain, and so we will investigate several different choices
to the less massive before this step). Finally, we place based on observations.
the star randomly on the surface of the sphere of radius For most of our simulations that include binaries, we
r . In this way we ensure that the clusters always match assume a binary system fraction, f = 0.5. We adopt
i bin
4 Farias, Tan & Chatterjee
Table 1
Parentclumpparameters
Σcloud Mcl(M(cid:12)) Rcl (pc) kρ φP,cl φP¯ φB σcl (km/s)
“Low-Σ”Clump 0.1 3000 1.159 1.5 2 1.31 2.8 1.71
“High-Σ”Clump 1 3000 0.367 1.5 2 1.31 2.8 3.04
Table 2
Initialconditionsforsimulationsets
Setname (cid:15) (cid:104)N∗(cid:105) fbin f(e) IMF IMS S.E. Commentinplots
equal mass 0.5 1500 0 – singlemass N N Singleequalmassparticles
single imf 0.5 4000 0 – Kroupa(2001) N N SinglestarswithIMF(NoSE)
binaries 50 0.5 4000 0.5 δ(e) Kroupa(2001) N N 50%binaries(NoSE)
fiducial 0.5 4000 0.5 δ(e) Kroupa(2001) N Y FiducialCase
binaries un 0.5 4000 0.5 uniform Kroupa(2001) N Y euniformdistribution
binaries th 0.5 4000 0.5 2e Kroupa(2001) N Y ethermaldistribution
segregated 0.5 4000 0.5 δ(e) Kroupa(2001) Y Y Masssegregated
binaries 100 0.5 4000 1 δ(e) Kroupa(2001) N Y 100%binaries
sfe 10 0.1 850 0.5 δ(e) Kroupa(2001) N Y SFE=10%
sfe 30 0.3 2500 0.5 δ(e) Kroupa(2001) N Y SFE=30%
sfe 80 0.8 6500 0.5 δ(e) Kroupa(2001) N Y SFE=80%
sfe 100 1.0 7300 0.5 δ(e) Kroupa(2001) N Y SFE=100%
Note. —Foreachofthesetsnamedinthefirstcolumn20simulationswereperformedforeachoftheclumpsparameterslistedinTable
1. SecondcolumnshowstheassumedSFE,thirdcolumnshowstheaveragenumberofstarspersimulation,fourthcolumntheprimordial
binary fraction, fifth column the eccentricity distribution function, column six shows the assumed IMF, column seven stands for whether
theclusterisinitiallymasssegregated(IMS).ColumneightshowsifStellarevolutionisincludedinthesetandlastcolumnisacomment
fromwhichwewillreferringtothesetinthegraphsforclarity.
a period distribution from the survey (Raghavan et al. We start using only single mass particles of m = 1M
i (cid:12)
2010) using a log-normal period distribution with mean with no primordial binaries and with a SFE of 50% in
of P = 293.3 years and standard deviation of σ = the set equal mass. Next, we include an IMF assum-
logP
2.28 (with P in days). We use a companion mass ratio ing the Kroupa (2001) distribution with a mass range
distribution (CMRD) of the form dN/dq ∝ q0.7, based from 0.01 to 100M(cid:12) defining the set single imf, again
onobservationsinyoungstarclusters(Reggiani&Meyer withnoinitialbinaries. Wetheninclude50%primordial
2011). The eccentricity distribution remains less well binaries with circular orbits and with other properties
constrained. While Duquennoy & Mayor (1991) found a described in §2.2, defining the set binaries 50. The
thermal distribution, i.e., f(e) = 2e for solar-type stars above three simulation sets do not include stellar evo-
in the solar neighborhood, a similar more recent study lution (SE). We define the fiducial simulation set by
(Raghavan et al. 2010) found a flat eccentricity distri- assuming a Kroupa IMF, SFE of 50%, with 50% of stars
bution for the same kind of stars. However, if binaries as primordial binary systems with initial circular orbits
formmainlyviadiskfragmentationwewouldexpectthat and with stellar evolution included.
theyarebornwithnearcircularorbits. Inordertomea- Next,wetestotherchoicesandparametersofthefidu-
sure how much binaries are affected by dynamical inter- cial model. We test two other eccentricity distributions,
actionsinthedifferentmodelsweadoptinitiallycircular a thermal eccentricity distribution, i.e., f(e)=2e in the
orbits for the eccentricities in our fiducial model. We setbinaries th,andauniformeccentricitydistribution
also investigate cases with initially thermal and uniform between0and1inthesetbinaries un. Anextremesce-
distributions of eccentricities. narioofmasssegregationistestedinthesetsegregated
inwhichstarsaresortedindescendingorderofindividual
2.3. Overview of the Cluster Models stellar mass from the center of the cluster. We also test
the extreme case in which all stars are binary systems
ForeachofthelowandhighΣclumps(seeTable1),we
(f =1) in the set binaries 100.
set up a stellar cluster as described above, i.e., assuming bin
We also carry out simulations with different SFE.
constant SFE(r) and a velocity dispersion profile equal
These simulations only differ from the fiducial set in
to that of the parent clump. Thus, the initial crossing
their SFE, i.e., the average number of stars per simu-
time (i.e., dynamical time) is defined by the properties
lation on each set increases with the SFE since we use
of the parent clump to be t ≡ R /σ , i.e., 0.663 and
cr cl cl the same parent clump of M = 3000M . The SFEs
0.118Myr for the Σ = 0.1 and 1gcm−2 cases, re- cl (cid:12)
cloud investigated are SFE = 10%, 30%, 80% and 100% and
spectively.
thesetsarenamedsfe 10,sfe 30,sfe 80andsfe 100,
Werun20realizationsforeachsetofinitialconditions,
respectively.
summarized in Table 2. The simulations are run for 20
Myr, varying only the random seed between simulations
2.4. Numerical Methods
in the same set, which affects the initial positions and
velocities of the stars as well as the IMF sampling. We follow the evolution of the star clusters for 20
We construct these sets of simulations starting from Myr utilizing the direct N-body integrator Nbody6++
thesimplestcasetotheonethatdefinesourfiducialcase. (Wang et al. 2015) which is a GPU/MPI optimized ver-
Star Cluster Formation from Turbulent Clumps 5
5 The kinetic energy of the stars is given by
3
T = M σ2, (13)
∗ 2 ∗
4
where σ is the one dimensional mass averaged velocity
dispersion. Assuming that the stars are born from the
gasfollowingthesamedispersionprofile,thenσisrelated
3
to that at the surface by Eq. 9.
Qi Replacing equations 9, 13 and 12 in eq. 10, and also
replacing R and σ by their expressions in eq. 4 and
2 φB=1 8 respectiveclyl, we obclt,asin
3(5−2k )(3−k ) 1
1 φB=2.φ8B=2 Qi = (8−3kρ)ρ2(kρ−1ρ)(cid:15)φB (14)
φB=5 0.51
0 φB=100 φB=20 Qi→ (cid:15) , (15)
0.0 0.2 0.4 0.6 0.8 1.0 where the arrow shows the relation using the fiducial
SFE
valuesfortheclump. ValuesofQ versusSFEareshown
i
Figure 1. The initial virial ratio Qi as a function of the SFE for different models in Figure 1.
fromtheclump. Solidlineshowstherelationforourfiducialvalue The dynamical state of the clusters also depends on
anddottedlinestherelationfordifferentvaluesofφB. the presence of magnetic fields in the initial clump, i.e.,
φ . In the absence of magnetic field support (φ = 1)
B B
the velocities needed for virial equilibrium are higher
sion of the classical and widely used direct integrator and stars formed from such gas will have higher val-
Nbody6 (Aarseth 2003). Nbody6++ has implemented ues of Q , e.g., even in the best case scenario with a
i
special regularizations to accurately follow the evolution SFE of a 100% we have Q (cid:39) 1.6 (and (cid:39) 3 for SFE of
i
of binaries and high order systems in the cluster being 50%). However, in the fiducial case with approximate
able to efficiently simulate star clusters with high binary equipartitionofenergydensityfromlargescalemagnetic
fractions with no loose of accuracy. For cases with stel- fields and turbulence (φ (cid:39) 2.8) a SFE of about 50%
B
larevolutionweusedtherecipeincludedinNbody6++ leadstoaclusterthatismarginallygravitationallybound
basedontheanalyticalmodelsforsingleandbinarystel- (Q (cid:39)1). Notethatthisvariationofφ alsocorresponds
i B
lar evolution developed by Hurley et al. (2000, 2002). to a variation in the virial parameter of the gas clump,
Thecodealsohasimplementedartificialvelocitykicksto α = 5(cid:104)σ2(cid:105)R/(GM), since for the fiducial case with
vir
emulate asymmetrical supernovae ejections. The magni- k =1.5 we have α =15/(4φ )→1.34 (see Appendix
ρ vir B
tude of the kicks are drawn from a Maxwell distribution A of MT03).
with σ = 265 km/s following the observations of Hobbs
et al. (2005) on pulsar proper motions. 3.2. The bound stellar cluster
As discussed in §3.1, star clusters born from turbulent
3. RESULTS clumps bounded by high pressure ambient environments
3.1. Initial dynamical state of the clusters can start with relatively high velocity dispersions. Their
virial ratio after gas expulsion will depend on the global
Before performing any simulation, from the assump-
SFE and the contribution of magnetic fields to the sup-
tionsdescribedin§2,wefirstderivetheinitialdynamical
portoftheparentclump,i.eφ . Therewillbesignificant
state of the clusters by characterizing their virial ratio, B
initial mass loss of the stars that are born unbound, oc-
i.e.,
curring on a timescale of a few crossing times. However,
T the gaseous clump is assumed to have a positive power
Qi=−Ω∗ (10) lawforthevelocitydispersionwithradius(seeEq. 9)and
so is more likely to be left with a central gravitationally
where T is the total kinetic energy of the stars and Ω is boundcore. Incontrast,arelaxedstarcluster(e.g.,with
∗
their total gravitational potential energy. A cluster with a Plummer profile) has a velocity dispersion profile that
Q < 1 is bound and Q = 0.5 is the value for a state decreases with radius. We thus expect differences in the
i i
of virial equilibrium. We assumed the gas was expelled early evolution of our clusters compared to those mod-
immediatelyafterthestarsformed, thusΩonlydepends eled with initial Plummer profiles by, e.g. Goodwin &
onthestarsinthecluster, i.e., Ω=Ω . Thepotentialof Bastian (2006); Baumgardt & Kroupa (2007); Pfalzner
∗
the stars is then & Kaczmarek (2013).
To measure the bound mass fraction, f , at each
Ω =−G(cid:90) Rcl(cid:20)M(r <Rcl)(cid:21)2dr timestep of cluster evolution. We construbcotuntdhe bound
∗ 2 r entity based on an accurate measure of the mean veloc-
0
ityoftheboundstarsandweselectallstarswithnegative
G(cid:90) ∞(cid:18)M (cid:19)2
− ∗ dr (11) total energy in the frame of reference of the bound clus-
2 r ter. The velocity of the bound cluster is not known a
Rcl
(cid:18) 3−k (cid:19)GM2 priori (although it is expected to be close to the zero ve-
=− ρ ∗. (12) locityofthereferenceframe),andthusthisissolvedinan
5−2k R
ρ cl
6 Farias, Tan & Chatterjee
1.0 presentedaretheaverageofthe20simulationsperformed
Measured at 2 t for each set. In each figure, the top four panels, (a) to
cr
(d), show the effects of gradually adding greater degrees
Measured at 20 Myr
of realism to make the fiducial model. Then the lower
0.8 cloud=0.1 gcm−2 four panels show the effects of different choices of initial
cloud=1 gcm−2 binary properties and degree of initial mass segregation.
These figures also show the evolution of the core radius
r (red dashed line), which is the density averaged dis-
c
0.6 tance from the density center of the cluster (see §15.2 in
Aarseth 2003).
d
un As expected, the clusters expand with time. The ex-
o
b pansionrateoftheouterLagrangianradiioftheclusters,
f
0.4 i.e.,oftheunboundstars,isdeterminedbytheinitialve-
locity dispersion of the parent clump, which is higher at
higher mass surface densities. Thus star clusters start-
ing from a clump with Σ =1gcm−2 are soon more
cloud
0.2 extended than the clusters forming from clumps with
Σ =0.1gcm−2 of the same mass and SFE, i.e., the
cloud
half-mass radius at 20 Myr of the first group is ∼ 20pc
compared to ∼10pc for the lower Σ case.
0.0 Initial expansion of the bound portion of the cluster
0.0 0.2 0.4 0.6 0.8 1.0
happens early within a few crossing times as the clus-
SFE ters relax to a virialized state. Then the later evolution
is affected by dynamical interactions between the stars
Figure 2. Bound mass fraction, fbound, as a function of SFE. (i.e., mass segregation, evolution of binaries and dynam-
Solid lines and points show results at early times (2tcr); dashed icalejectionofstarsfromunstablemultiplesystems)and
lines and open points show results at the end of the simulations
(20 Myr). The values are the mean averages for the simulation mass loss resulting from stellar evolution. The relative
sets (each of 20 realizations) with cases shown for Σcloud = 0.1 importance of these effects can be gauged by examining
(green)and1gcm−2(red). Errorbarsshowthevaluesbetweenthe the sequence of panels from (a) to (d) in Figures 3 and
25thand75thpercentiles. Dynamicalevolutionleadstoageneral
4. Thelaterstageexpansionoftheboundclusterisneg-
decreaseoff duringthefirst20Myrsoftheseclusters.
bound
ligible in the case of equal mass stars. The model with
iterative calculation. We start by selecting all stars with
anIMFbutonlysinglestarsundergoesmasssegregation
negative energy inside the half mass radius of the full
that leads to noticeable expansion after about 6 Myr in
cluster and then iterate until the members between iter- the case ofΣ =0.1gcm−2 and afterabout 1 Myrin
cloud
ations do not change by more than two members. This the case of Σ =1.0gcm−2.
cloud
method,called“snowballing”,isdescribedinSmithetal.
Note that before adding binaries and stellar evolution
(2013b).
in our models, evolution of the star clusters with high
Figure 2 shows the bound fractions measured at 2
Σ wouldbeexactlythesameasthosewithlowΣ
cloud cloud
crossing times and after 20 Myr for clusters with dif-
afterproperlyscalingfortheinitialsizeandcrossingtime
ferent SFE for simulations with a parent clump with
(see Aarseth & Heggie 1998). However, the characteris-
Σ = 0.1 and 1.0gcm−2. The initially positive ra-
cloud tic timescales introduced by binaries (e.g., at their typi-
dial gradient of velocity dispersion of the star clusters
calorbitalseparation)andbystellarevolutionbreakthis
causes stars in the outer parts to leave first, while the
self-similarity.
central core can remain bound. The mass fraction of
When binaries are added, their presence leads to an-
this remnant bound core depends on the initial global
other potential source of expansion, since their binding
SFE. Later dynamical evolution and internal mass loss
energy couples with the internal energy of the stellar
of its members due to stellar evolution leads to a slower
cluster, leading to a change of kinetic energy in each
decrease in mass of the bound core over time.
interaction (Heggie 1975; Hills 1975). However little dif-
TheresultsshowninFigure2arethoseforthefiducial
ference appears when moving from single stars to 50%
case, i.e., with φ =2.8. As discussed in §3.1, the initial
B binaries, even in simulations with the high Σ initial con-
virialratioofthestarclustersdependssensitivelyonthis
dition that can suffer more interactions. As we will see
value: ahighervalueofφ shiftsthetrendshowninthis
B in §3.4, the initial densities of these models are not high
figure upwards so that even clusters with SFE = 10%
enough and/or do not last long enough to give binaries,
may retain a significant bound core if φ is sufficiently
B onaverage,thechancetointeractsignificantlywithother
high.
stars.
The inclusion of stellar evolution causes the cluster to
3.3. Global evolution
expandevenmore. Stellarevolutionstartsbecomingim-
Here we explore the evolution of the structure of the portantafterafewMyr,whenthefirstmassivestarslose
clusters with time. The fiducial case has SFE of 50% mass by stellar winds and then explode as supernovae.
and Q = 1.02, slightly above the criterion for global The supernova explosions may cause stars to be ejected
i
boundedness. Therefore, initial expansion and some ini- (e.g., as fast runaway stars) either by the destruction of
tial mass loss is expected. We show the evolution of tightbinaries, velocitykickscausedbytheasymmetrical
the Lagrangian radii with time in Figures 3 and 4 for explosion, or both combined effects (see the outer two
Σ =0.1and1gcm−2,respectively,wherethevalues lines in Figures 3 and 4). We focus on the ejection of
cloud
Star Cluster Formation from Turbulent Clumps 7
102
(a) Single equal mass particles (b) Single stars with IMF (No SE)
101
)
pc 100
(
s
s
a
m10-1
al
t
o
% t10-2
0 102
9 (c) 50% binaries (No SE) (d) Fiducial Case
%,
80 101
%,
0
%, 7 100
0
6
%, 10-1
0
5
%, 10-2
40 102
%, (e) e uniform distribution (f) e thermal distribution
0
%, 3 101
0
2
%, 100
0
1
e
h10-1
t
r
o
f
dii 10-2
ra 102
n (g) Mass segregated (h) 100% binaries
a
gi
an 101
r
g
a
L
100
10-1
10-2
10-1 100 101 10-1 100 101
time (Myr) time (Myr)
Figure 3. Average Lagrangian radii evolution for different set of simulations for star clusters born from a parent clump with Σ =
cloud
0.1gcm−2. We show the Lagrangian radii for the 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80% and 90% masses (black solid lines). Red
dashed lines are the core radii defined in Aarseth (2003). Gray shaded areas represent the regions below the 50%, 95% and 100% mass
radiusoftheboundcluster.
runaway stars below in §3.6. For now we see that clus- effects.
ter expansion is increased by this effect and also by the In Figure 5 we compare the evolution of several pa-
fact that the potential well of the cluster is made more rameters of the different sets of simulations for clus-
shallow with the loss of mass through stellar winds and ters born with different initial mass surface densities:
supernovae. However, the loss of runaway stars does not Σ =0.1gcm−2 ontheleft; Σ =1gcm−2 onthe
cloud cloud
affect the global evolution of the cluster too much, even right. In the first row we show evolution of the bound
in the most extreme case with f = 1. Finally, the mass fraction f . Only for the purposes of this fig-
bin bound
lower four panels in these figures show that variations of ure, to show the timescale on which initially unbound
binary orbital properties, degree of initial mass segrega- stars leave the cluster, we count all stars inside the 95%
tion or primordial binary fractions have relatively minor radius of the bound cluster as also being bound. Here,
8 Farias, Tan & Chatterjee
102
(a) Single equal mass particles (b) Single stars with IMF (No SE)
101
)
pc 100
(
s
s
a
m10-1
al
t
o
% t10-2
0 102
9 (c) 50% binaries (No SE) (d) Fiducial Case
%,
80 101
%,
0
%, 7 100
0
6
%, 10-1
0
5
%, 10-2
40 102
%, (e) e uniform distribution (f) e thermal distribution
0
%, 3 101
0
2
%, 100
0
1
e
h10-1
t
r
o
f
dii 10-2
ra 102
n (g) Mass segregated (h) 100% binaries
a
gi
an 101
r
g
a
L
100
10-1
10-2
10-1 100 101 10-1 100 101
time (Myr) time (Myr)
Figure 4. SameasFigure3forstarclustersbornfromaclumpwithΣ =1gcm−2.
cloud
and in all panels in this figure, values are the medians of tion.
eachsetofsimulationswithparametersgiveninTable2. It takes about 1.5 t for initially unbound stars to
cr
Also shown in the top row is a shaded area representing leave the bound cluster, leading to f decreasing to
bound
the loss of mass due to stellar evolution for all stars in about0.7. Thevelocityprofileoftheclustershaveapos-
the fiducial simulations (including unbound stars). The itiveslope,i.e.,higherspeedsintheoutskirts,thiscauses
second row shows the evolution of the core radii, r , and outer stars to be more likely to be unbound, with prac-
c
the half-mass radii, r . The third row shows the evolu- tically no chance of interacting with others. These stars
h
tion of the effective number density, n , i.e., the number leave the cluster with a velocity dispersion determined
s
ofsystems(abinaryiscountedasonesystem)insidethe by the parent clump, i.e. σ . We refer to these as un-
cl,s
volume defined by r . The fourth row shows the evolu- bound stars, distinguishing from the ejected stars that
h
tion of the velocity dispersion measured inside r , while escape later due to dynamical ejections. After the first
h
thefifthrowshowstheevolutionofthetotalbinaryfrac- ∼1.5t ,allinitiallyunboundstarsleavetheclusterand
cr
Star Cluster Formation from Turbulent Clumps 9
Single equal mass particles 50% binaries (No SE) e uniform distribution Mass segregated
Single stars with IMF (No SE) Fiducial Case e thermal distribution 100% binaries
1.0
0.9 Σ =0.1 g/cm2 Σ =1 g/cm2
cloud cloud
0.8
d
n
ou0.7
b
f
0.6
0.5
0.4
101
]
c
p 100
[

h
r

,
rc10-1
10-2
104
]
3c 103
p
/
s
r 102
a
t
s
[
s101
n
100
4.0
3.5
3.0
]
s 2.5
/
m
2.0
k
[
1.5
σ
1.0
0.5
01..00
0.9
0.8
n0.7
bi
f
0.6
0.5
0.4
10-1 100 101 10-1 100 101
t [Myr] t [Myr]
Figure 5. TimeevolutionofpropertiesoffiducialclustersformingwithSFE=50%. LeftcolumnshowsclustersformingfromaΣ =0.1
cloud
gcm−2environment;rightcolumnfromaΣ =1gcm−2environment. Thelinesineachpanelshowmedianvaluescalculatedfromthe
cloud
20simulationsperformedforeachset. Toprowshowsthefractionofboundmassintheclusterrelativetotheinitialmass,whereinthis
figureunboundstarsinsidethe95%massradiusoftheboundclusterarekepttoshowthetimescaleoftheirescape. Herewealsoshowthe
fractionoftotalstellarmassinthefiducialsimulationsthatremainsafteraccountingforstellarevolutionmassloss: grayshadedregion
showstheareabetween25thand75thpercentiles. Secondrowshowstheevolutionofcoreradius(rc)andhalfmassradiusrh forallthe
starsinthesimulation. Thirdrowshowstheaveragenumberdensityofsystems(ns),wherepersystemswerefeertosinglesandbinaries,
measuredinsidethevolumedefinedbyr . Fourthrowshowtheevolutionofthevelocitydispersionmeasuredinsider , andbottomrow
h h
showstheevolutionoftheglobalbinaryfraction.
10 Farias, Tan & Chatterjee
later evolution is determined by dynamical interactions 103
and stellar evolution. All stars
Simulations with equal mass stars essentially do not Bound stars
lose further members. With an IMF, then mass segrega- Σ =0.1g/cm2
cloud
tion does lead to some additional mass loss. When in- 102 Σcloud=1g/cm2
cluding 50% primordial binaries mass loss at later times
ismoderatelyenhanced. Addinginstellarevolution,i.e., h,i
r
in the fiducial model, continues this trend, with a final /
value of fbound (cid:39) 0.5. These trends are mirrored in the rh,f
expansion of the clusters. Variations of binary orbital 101
propertiesorprimordialbinaryfractionsareseentohave
relatively minor effects.
Models with full initial mass segregation show some
differences. In the case with Σ = 0.1gcm−2, the
cloud 100
numberdensitiesatthecenterarequitelowinitiallyand 102 All stars
the core of the cluster contracts significantly. After this
Bound stars
contraction, the number density is raised in the core,
Σ =0.1g/cm2
which then later expands quite rapidly. Even though cloud
Σ =1g/cm2
number densities of these clusters are never too high, cloud
the few interactions that do occur are enough to expand
the cluster and the evolution of the 50% mass radius is
101
determined by these interactions.
)
Ingeneral,thestarclusterspresentedhereexpandcon- c
p
siderably regardless of the different parameters of the (
simulations. The amount of expansion depends on the h,f
r
SFE and the initial cluster density. The top panel of
Figure6showstheratiobetweenthehalf-massradiusat
theendofthesimulations,r ,andthehalf-massradius 100
h,f 0.0 0.2 0.4 0.6 0.8 1.0
atthestart,rh,i. Thedifferencebetweenthemodelswith SFE
highandlowinitialdensitiesareexplainedmainlybythe
Figure 6. Thesizeofstarclustersat20Myrasafunctionofthe
initialvelocitydispersionoftheparentclump. Starsthat
SFE.Toppannelshowsthefinalhalf-massradii,r ,comparedto
arebornunboundinthedensercaseescapewithahigher theinitialhalf-massradii,r . Bottompanelshowhs,rf inphysical
h,i h,f
typical velocity than the low dense case, causing the dif- units. Filled circles show measurements using all the stars; open
ferences in the expansion. However, when considering circlesusingonlytheboundcluster. Valuesaremediansoverthe20
simulationsperformedforeachsetanderrorbarsshowstheregion
theboundpartofthecluster,theactualsizeofstarclus-
betweenthe25thand75thpercentiles.
ters at 20 Myr is similar regardless the initial density, as
shown in the bottom panel of Figure 6. Differences be-
We now estimate the perturbation encounter rate of
tweenthesizesoftheboundclustersarisewhentheSFEs
a binary of a given semimajor axis a in our model star
are low. This is due to the fact that the crossing times
clusters. We first derive the rate assuming stars move
oftheseboundsystemsarelarge(t ≈30Myr)andthey
cr without significant deflection, then include the effects of
have not yet achieved an equilibrium distribution by 20
gravitational focusing. If we assume that the cluster has
Myr. Thus low SFE clusters are still in the first phase
onlysinglestarsandbinaries,themeanmasspersystem
of their expansion. Regardless of the initial density, the
is (cid:104)m (cid:105) = (1 + f )(cid:104)m (cid:105). If there exist higher order
s bin i
final size of the bound systems depends mainly on the
multiples, then (cid:104)m (cid:105) will be higher, however these are
s
initial SFE: low SFE results in a more extended bound
not included as initial conditions in our models and we
system.
will see that interactions are typically at a relatively low
rate so that such multiples do not form in significant
3.4. The effects and evolution of binaries
numbers during the dynamical evolution of the clusters.
Our modeling includes a full treatment of binaries, so The mean rate of interactions that are able to modify
we are able to examine their effects and evolution in de- thepropertiesofabinary,Γ ,isproportionaltothecross
b
tail. A binary will be significantly perturbed by an ex- section defined by b, i.e., πb2, multipled by the number
ternal star (or multiple) if the potential energy of the density of perturbing systems n and a typical velocity
s
encounter is similar to that of the initial binary, i.e., in the cluster, i.e., the 1D velocity dispersion, σ. Thus
the interaction rate for binaries with a given semi-major
Gm m G(m +m )m
E =− 1 2 ∼− 1 2 s, (16) axis a is
bin 2a b
(cid:18)(cid:104)m (cid:105) (cid:19)2
where b is the closest approach of the perturber of mass Γ =4π s a n σ. (18)
m , and m and m are the primary and secondary b µ s
s 1 2
masses of the binary, respectively. Therefore, defining
As we show in Figure 5 the effective number density
µ≡m m /(m +m )asthereducedmassofthebinary,
1 2 1 2 in our model clusters quickly falls from initial values of
the closest approach needed to affect binary properties
∼ 103 to 104 stars/pc3 (depending on the initial en-
is:
vironment mass surface density) to values similar to 1
m
b∼2a s. (17) stars/pc3 at 20 Myr, in our fiducial case. The typical
µ