The fragmentation of expanding shells III: Oligarchic accretion and the mass spectrum of fragments
James E. Dale, Richard Wunsch, Rowan J. Smith, Anthony Whitworth, Jan Palous
aa r X i v : . [ a s t r o - ph . GA ] O c t Mon. Not. R. Astron. Soc. , 1–12 (2006) Printed 9 November 2018 (MN L A TEX style file v2.2)
The fragmentation of expanding shells III: Oligarchicaccretion and the mass spectrum of fragments
James E. Dale ⋆ ,Richard W¨unsch , ,Rowan J. Smith ,Anthony Whitworth ,Jan Palouˇs Astronomical Institute, Academy of Sciences of the Czech Republic, Bocni II 1401/2a, 141 31 Praha 4 School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff, CF24 3AA Institut f¨ur Theoretishce Astrophysik, Universit¨at Heidelberg, Albert– ¨Uberle–Str. 2, Heidelberg, Germany
ABSTRACT
We use SPH simulations to investigate the gravitational fragmentation of expandingshells through the linear and non–linear regimes. The results are analysed using spher-ical harmonic decomposition to capture the initiation of structure during the linearregime; the potential-based method of Smith et al. (2009) to follow the developmentof clumps in the mildly non-linear regime; and sink particles to capture the propertiesof the final bound objects during the highly non-linear regime. In the early, mildlynon–linear phase of fragmentation, we find that the clump mass function still agreesquite well with the mass function predicted by the analytic model. However, the sinkmass function is quite different, in the sense of being skewed towards high-mass ob-jects. This is because, once the growth of a condensation becomes non-linear, it tendsto be growing non-competitively from its own essentially separate reservoir; we callthis Oligarchic Accretion.
Key words: stars: formation, ISM: HII regions
Bubbles and the shells they sweep up are ubiquitousfeatures of the interstellar medium (ISM), on a widerange of mass and size scales. Churchwell et al. (2006) andChurchwell et al. (2007) observed over 600 infrared shellswith Spitzer on size scales of 0 . . ∼ ⋆ E-mail: [email protected] (JED) consensus on the mass function of stars produced by thisprocess. Whitworth et al. (1994b) assume that the most un-stable mode in the shell can be used to characterize the massfunction, concluding that shell fragmentation should prefer-entially produce high-mass stars. If this is correct, it stronglysupports the self-propagating star formation model, sinceeach shell produces more than enough O-type stars to trig-ger the next stellar generation. However, W¨unsch & Palouˇs(2001) perform a more detailed analysis under the assump-tion that all unstable wavelengths are represented in themass function. They derive a fragment mass function fromthe dispersion relation for the growth rate of perturbationsof different angular wavenumbers, and conclude that shellfragmentation should produce a power–law mass functionnot very different from the canonical Salpeter function.Unambiguous observational evidence for shell fragmen-tation is difficult to obtain. Churchwell et al. (2007) findthat ∼
18% of their sample of 269 bubbles in the innerMilky Way show evidence of triggered star formation, butattribute this to the shells over–running pre–existing densityanomalies, rather than to fragmentation of the shells them-sleves. Observations of shells driven by HII regions by, forexample, Zavagno et al. (2006), Deharveng et al. (2006) andDeharveng et al. (2008) often detect populations of youngstars in the peripheries of bubbles and the stars are of-ten massive. Herschel observations of Sh2–104 (Rod´on et al.2010) and RCW120 (Zavagno et al. 2010) seem to imply c (cid:13) James E. Dale, Richard W¨unsch, Rowan J. Smith, Anthony Whitworth, Jan Palouˇs that shell fragmentation produces small numbers of largefragments. However, observations by Dawson et al. (2008)of the mass function of CO clumps in the Carina Flare su-pershell reveal a power–law mass spectrum. It is thus notclear from either an observational or a theoretical point ofview how self–gravitating shells should fragment.The present paper is one of a series aimed at resolvingthese issues by performing detailed numerical simulationsof fragmenting shells. In the first paper (Dale et al. (2009);hereafter Paper I) we studied the fragmentation of expand-ing shells in the linear regime using a Smoothed ParticleHydrodynamics (SPH) code and an Adaptive Mesh Refine-ment (AMR) code. We studied the growth of perturbationsin the shell surface density and compared the results of thetwo codes with each other and with the thin–shell model.We found that the agreement between the two radically dif-ferent numerical schemes was excellent. However, we alsofound that the boundary conditions applied to the shell wereof crucial importance in determining the mass-spectrum ofperturbations. A shell expanding in a vacuum gets thickeras it expands, and this thickening suppresses fragmentationat short wavelengths, relative to the predictions of the thinshell model. Conversely, if the thickening of the shell is pre-vented by a confining pressure, the shorter wavelengths be-come unstable too and, for a certain finite value of the pres-sure, fragmentation proceeds in good agreement with thethin shell model.The second paper in this series, Wunsch et al. (2010)and hereafter Paper II, presents an alternative to the thin–shell model which allows the boundary conditions on theshell inner and outer surfaces to be included and thus theeffects of external pressure to be treated analytically. Wefind that pressures higher than that required to keep theshell thickness constant extend the range of unstable wave-lengths towards lower values. This should in principle resultin a fragment mass function with a steeper slope and weinvestigate this possibility below.In this paper, we extend some of the simulations of Pa-per I into the non-linear regime and allow the fragmentationof the shells to proceed to the formation of bound com-pact objects. These can be thought of either as protostarsor protoclusters, depending on their masses, allowing us (i)to derive a mass function for the stellar/cluster populationsformed by fragmenting shells; and (ii) to evaluate the abilityof the thin–shell model to predict these mass functions. InSection 2, we discuss the thin–shell model and the meansby which mass functions may be inferred from it. In Section3, we explain our numerical methods and in Section 4, webriefly outline the simulations conducted. Our results arepresented in Section 5 and discussed in Section 6. Finally,we summarize our conclusions in Section 7.
For a self–gravitating shell with instantaneous radius, ex-pansion velocity and mean surface density given by R , V andΣ , the infinitesimally thin shell model yields an expressionfor the growth rate ω of perturbations in the surface density of spherical harmonic wavenumber l , ω ( l ) = − AVR + r BV R − c s l R + 2 πG Σ lR , (1)where c s is the sound speed inside the shell, and A and B areconstants whose values depend on whether the shell sweepsup mass as it expands. In the case of a shell expanding into avacuum, A = , B = . In the shells studied in Paper I andin this paper, the expansion velocity V of the shell is smallcompared with c s and so the first two terms in Equation 1are small.In Wunsch et al. (2010), we derived an alternative ex-pression for the dispersion relation which drops the assump-tion that the shell is thin and takes into account the effectof the external pressure, P EXT ω ǫ = − V ǫR + ((cid:18) V ǫR (cid:19) + 3 πG Σ ǫ " r cos − ( z /r )( r − z ) / − z ( r − z ) + 10 P EXT c S l π ǫR (2 P EXT + πG Σ ) − c S l π ǫR (cid:27) / , (2)where z and r are respectively the vertical half–thicknessand azimuthal radius of a perturbation and ǫ is a normalisa-tion parameter, described in Wunsch et al. (2010). In orderto be useful in the study of star formation, these dispersionrelations must be converted into fragment mass functions.Beginning from the thin–shell model, Whitworth et al.(1994b) neglect the shell expansion terms and consider onlythe interplay between gravity and gas pressure within theshell and the masses of fragments formed. They considerperturbations of radius r , compute the acceleration g of ma-terial inside them, and write the timescale on which such afragment condenses out of the shell as t frag ≈ (cid:18) rg (cid:19) / ≈ (cid:20) G Σ r − (cid:16) c s r (cid:17) (cid:21) − / . (3)Whitworth et al. (1994b) then find the value of r that min-imises Equation 3 and take the mass corresponding tothis fragment size as representative of the mass function.Dale et al. (2007) simulate the early stages of the fragmen-tation of a shell driven by an expanding HII region andfind reasonable agreement with Whitworth et al. (1994b) interms of the time and shell radius at which fragmentationbegins and the mean masses of fragments formed.W¨unsch & Palouˇs (2001) introduce a more sophisti-cated analysis in an attempt to derive the mass spectrumof fragments, but their derivation is not correct, since thereis an error in their Equation 52. Although we import someof their results, we provide here a different semi–analyticderivation.We seek an expression for the number of fragments dN existing at a time t in the range of masses [ m f , m f + d m f ]and we make the assumption that the right–hand side of thisexpression can be separated into two parts; the fragmenta-tion integral (now written as a function of mass) I f ( m f , t )),which expresses how much perturbations of a given masshave grown at time t , and a geometrical function G ( m f )d m f c (cid:13) , 1–12 he fragmentation of expanding shells III: Oligarchic accretion and the mass spectrum of fragments Figure 1.
Result of an experiment of covering a square withcircles uniformly drawn in mass from (cid:2) . × − , . × − (cid:3) up to a covering factor of 0.8. which encodes how many fragments of a given mass will fiton the shell. We may then write dN ( m f , t ) = I f ( m f , t ) G ( m f ) d m f . (4)The fragmentation integral is simply defined as I f ( m f , t ) = Z tt ω ( m f , t ′ )d t ′ , (5)where ω ( m f , t ′ ) is the instantaneous growth rate of the modecorresponding to m f and t is the time at which fragmenta-tion begins.The derivation of the function G ( m f )d m f is the point atwhich we depart from the methodology of W¨unsch & Palouˇs(2001). It is not clear to us that an analytic expression for G ( m f )d m f may be found since, on the reasonable assump-tion that fragments may not overlap each other, the numberof objects of any given wavenumber that may be accommo-dated by the shell depends on how much space has alreadybeen consumed by other fragments. We therefore adopteda Monte Carlo approach in an attempt to derive a semi–analytic approximation to G ( m f )d m f .We assume that fragments are circles. We took a unitsquare (the shape of the area to be filled is largely irrelevant,since circles do not tesselate) and populated it with circleschosen uniformly in [ m min , m max ], m f being the mass of acircle and proportional to its area. We forbade overlappingand continued until the square was covered up to a factor f by circles. We then constructed histograms of the popu-lation of circles. In Figure 1 we show the result of one suchexperiment in which f = 0 . m min , m max are 1.26 × − and 7.9 × − respectively. In Figure 2, we show the powerlaws resulting from hundreds of such experiments in whichwe varied f . Varying the covering factor affects the low–massend of the mass spectrum, which turns over at small massesdue to underrepresentation of small circles. Achieving highercovering factors requires more small circles and pushes theturnover towards smaller masses. At the high–mass end ofthe mass function, the mass function becomes noisy. How-ever, in between these limits, the slope of the function isvery robust. We find that the population of circles can bewell approximated by d N ∝ m − . d m f .Armed with this result we may write -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5-5 -4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 l og d N / d m log m Figure 2.
Mass functions of circles, with covering factors f of0.70 (plus signs), 0.75 (crosses) and 0.80 (triangles). The straightline has a logarithmic slope of -2.0. dN ( m f , t ) = I f ( m f , t ) m − . d m f . (6)An analytic expression for the fragmentation integralitself can be derived for the infinitesimally thin shell and isgiven in Appendix A in a dimensionless form that can be ap-plied to any thin shell. However, it is somewhat cumbersomeand we compute fragmentation integrals from numericallyintegrating Equation 5 using analytic wavenumber growthrates for the thin–shell and PAGI models and compare theresult to that computed from our SPH simulations. We make use of the same SPH code used in Paper I, a variantof that described in Benz (1990) but more recently updatedand described in Bate et al. (1995). The fluid equationsare solved using the SPH technique implementing the stan-dard artificial viscosity prescription described in Monaghan(1992), with α = 1, β = 2. The self gravity of the gas isincluded using a binary tree. Crucially for these calculation,the code allows very high density regions to be replaced bypoint–mass sink particles, as described in Bate et al. (1995).Once the density of a particle exceeds a threshold, setto 10 − g cm − in these simulations, it and its ≈
50 neigh-bours are considered candidates for sink creation. In orderfor a dense particle and its neighbours to be replaced by asink, the group of particles must (a) constitute more than athermal Jeans mass; (b) be contracting; (c) be bound. Gasparticles may subsequently be accreted by sink particles ifthey pass within the sink particle’s accretion radius and (a)are bound to the sink; (b) are more bound to that sink thanto any other sink; (c) do not have enough angular momen-tum to achieve a circular orbit around the sink. The massresolutions in the medium– and high–pressure simulationspresented here are 1.67M ⊙ and 0.4M ⊙ respectively. We usea sink accretion radius of 0.05pc which ensures that the ac-cretion flows at the accretion radius are always supersonic,so that accretion of gas particles does not produce spuriousdrops in pressure. We repeated part of the medium–pressurecalculation using an accretion radius of 0.025pc, but foundthat this made no difference to the results. The code alsoallows for corrections to be made to the pressure and den- c (cid:13) , 1–12 James E. Dale, Richard W¨unsch, Rowan J. Smith, Anthony Whitworth, Jan Palouˇs sity gradients in the neighbourhood of sink particles to com-pensate for the fact that there are few or no gas particlesinside the accretion radius, but we found that they had noeffect on our results. In these calculations, the gas flows arenon–turbulent and thus rather quiet and smooth, and hencecomparatively easy to model.To achieve greater accuracy, the sink particles are notincluded in the binary tree and all gravitational forces in-volving them are computed by direct summation. To preventsink particles in very close proximity to each other acquiringvery short timesteps and stalling the simulation, we allowedpoint masses to merge if they passed within 0.05pc of eachother and were mutually bound. During the course of thesimulations, there are a few tens of such mergers, too few tostrongly influence the mass functions produced.We use the same spherical harmonic decompositiontechnique discussed in Paper I to analyse the growth of per-turbations in the shell surface density, which we compare tothe thin shell model. The thin–shell model uses linearizedequations and is therefore only guaranteed to be valid inthe linear phase of fragmentation, where surface density per-turbations are less than or comparable to the mean surfacedensity. In this work we are interested in the transition of theshell into the non–linear regime and the contraction of suchperturbations into bound objects. Apart from this phase offragmentation being non–linear, there also comes a pointfor each contracting fragment when its wavelength becomescomparable to the thickness of the shell and fragmentationtransitions from a two–dimensional to a three–dimensionalprocess. In addition, the wavelengths and wavenumbers offragments changes as they collapse, and the surface densityassociated with the fragment grows, so that the associationof a given wavenumber with a given mass may become in-valid. For these reasons, it is not clear that the thin–shellor PAGI analysis may be extended into this regime of theshell’s evolution and we therefore turn to a numerical study.We use two different techniques to follow the fragmen-tation process beyond the linear phase. Once fragmentationhas become strongly non–linear and bound objects are col-lapsing on their local freefall times, we make use of SPHsink particles (Bate et al. 1995). However, in the regime be-tween the linear and strongly non–linear phases of the shellevolution, before any fragments have become gravitationallybound, we require a robust clump–finding method.Three dimensional clump–finding in complex densityfields is a difficult task and several numerical methodshave been developed to facilitate it, e.g. the CLUMPFIND(Williams et al. 1994) algorithm. In a particle–based codesuch as SPH, progress can be made in identifying perturba-tions in the density field by identifying contiguous regionsof fluid whose density exceeds some threshold. Constructinga mass function is problematic however, since the chosendensity threshold effectively selects the masses of the frag-ments. This problem can be solved by insisting that frag-ments must be bound. By selecting the densest particles asseeds for clumps and continually adding adjacent particlesuntil the total energy of the composite object becomes posi-tive, it is possible to identify a population of bound objects.This method was used by Dale et al. (2007), but is oflittle use here, since it is of no help in the mildly non–linearregime. We find that a given object goes very rapidly frombecoming bound to forming a sink particle, so that this I f I f Figure 3.
Comparison of the PAGI low–pressure (1 × − dynecm − , green), PAGI medium–pressure (1 × − dyne cm − ,black), PAGI high–pressure (5 × − dyne cm − , magenta),thin–shell (blue) and numerical (red) medium–pressure fragmen-tation integrals computed up to a time of 8.80 Myr. method reveals little that cannot be learned from the sinkparticles themselves. We therefore make use of the clump–finding method of Smith et al. (2009) which identifies struc-ture using peaks in the gravitational potential field, insteadof in density. There are two advantages to such a method.First, the gravitational potential distribution in the cloud isconsiderably smoother than the density distribution, sincedensity fluctuations that do not carry sufficient mass can-not contribute significantly to the potential field. Second, thestrength of the gravitational potential determines whether aclump will collapse and how mass will flow. The structuresidentified by this algorithm are called p-cores to distinguishthem from conventional density cores.The potential–clumpfinding algorithm works in a simi-lar manner to CLUMPFIND and can be applied directly tothe SPH data. The SPH particle with the deepest gravita-tional potential forms the head of a clump, then the particlewith the next deepest potential is either assigned to thesame clump if it is an SPH neighbour of the head particle,or forms a new clump if not, and so on. Clumps are defineddown to either a minimum positive potential, or the low-est contour which they share with a neighbouring clump.Unlike the traditional CLUMPFIND algorithm, we use con-tour levels primarily to define the level at which potentialclumps join, rather than to distinguish clumps from noise,and so our contours are numerous and finely spaced. Thishas the effect of subtracting the background potential. Thisis necessary as gravity is a long range force, and so is af-fected by both the mass inside the p–core and surroundingit; hence we must remove the background to obtain the neteffect on the mass within. P–cores, therefore, represent thelocal maximum above the surrounding background. c (cid:13) , 1–12 he fragmentation of expanding shells III: Oligarchic accretion and the mass spectrum of fragments − d N ( M ) / d M d N ( M ) / d M ⊙ )M (M ⊙ ) Figure 4.
Comparison of the PAGI low–pressure (1 × − dyne cm − , green), PAGI medium–pressure (1 × − dynecm − , black), PAGI high–pressure (5 × − dyne cm − ,magenta), thin–shell medium–pressure (blue) and numericalmedium–pressure (red) mass functions, derived from Equation6, computed up to a time of 8.80 Myr. The dashed line has alogarithmic slope of -2.25. . . . N o . o b j ec t s ⊙ ) Figure 5.
Comparison of clump mass functions generated byimposing an artificial surface density profile computed from thefragmentation integral of a given SPH dump on a uniform shell(red) with that computed from analysing the same dump directly(blue).
In Paper I, we presented the results of simulations ofexpanding momentum–driven shells with two differentboundary conditions. The shells had a mass of 2 × M ⊙ ,an initial radius of 10 pc and an initial velocity of 2.1kms − , such that they expand to a maximum radius of ≈ . × − dyne cm − to the inner and outer faces ofthe shell such that its thickness remained approximatelyconstant and found that this resulted in much betteragreement with the thin shell model. However, in Paper I,we followed the shell evolution only in the linear regime inwhich the thin–shell model is most likely to be applicable.In this paper, we allow the simulations of Paper I to runfurther, into the non–linear regime, using only the SPHsimulations and not the AMR simulations, since both sinkformation and clumpfinding are much more difficult inthe latter. We find that the fragmentation process in thenon–pressure confined simulation from Paper I proceedstoo slowly to form a significant number of sink particlesbefore the shell has contracted back to its original size andwe do not examine it further.In Paper II, we considered the effects of higher pres-sures and found that they result in additional shorterwavelengths becoming unstable, which should producemore low–mass fragments. In this work, we concentrateon the pressure–confined model from Paper I (which werefer to as the medium–pressure model), and we also studythe mass function produced by the same shell, but with aconfining pressure of 5 . × − dyne cm − , which we alsostudied in Paper II and will refer to as the high–pressuremodel. In Figure 3 we plot the fragmentation integrals computedfor the medium–pressure run from the thin–shell model,for three different pressures using the PAGI model, andthat derived from the spherical–harmonic analysis of themedium–pressure SPH simulation, all computed at a time of8.80Myr. The analytic fragmentation integrals are obtainedby integrating Equation 5, using Equation 1 to compute ω in the thin–shell case, and Equation 2 to compute ω ǫ in the PAGI case. The numerical fragmentation integralclearly agrees much better with the relevant PAGI modelthan with the thin–shell model. Importing the resultsfrom Section 2, we then compute mass functions from thefragmentation integrals, shown in Figure 4. We again seethat the mass function computed from the PAGI model isin much better agreement with that computed from theSPH medium–pressure calculation than is the thin–shellmass function. We find that, at high masses, the theoreticalmass functions are power laws with slopes of -2.25, in goodagreement with the canonical Salpeter slope of -2.35.Before proceeding with clumpfinding, we sought to es-tablish a connection between the two–dimensional spectralanalysis technique and the intrinsically three–dimensionalconcept of a clump defined by its gravitational potential.We selected an output file from the early stages of themedium–pressure simulation and computed the correspond-ing fragmentation integral. We then used the integral to c (cid:13) , 1–12 James E. Dale, Richard W¨unsch, Rowan J. Smith, Anthony Whitworth, Jan Palouˇs . . . . . Σ ( g c m − ) Σ ( g c m − ) Figure 6.
Comparison of mean (solid line) and maximum per-turbed (dashed line) surface density in the medium–pressure cal-culation. construct a perturbed surface density profile which weimposed on a uniform shell by iteratively adjusting theSPH particle masses in a manner similar to that describedin Paper I. Finally, we computed a clump mass functionfrom the artificial shell and compared it to that computedfrom the original output file from the medium–pressuresimulation. As shown in Figure 5, the mass functionsgenerated in this way are in reasonable agreement, whichimplies that the spectral analysis technique is a legitimateway of describing shell fragmentation, at least in the earlylinear regime.
In Paper I, we performed tests in which we insertedmonochromatic density perturbations into the shell andcompared their evolution with that predicted by the thinshell model. We found that the agreement was very gooduntil the point when the perturbation in the surface densitybecame equal to the mean surface density, after which theperturbation growth became non–linear, departing stronglyfrom the model. In Figure 6, we plot the mean shell sur-face density (solid line) and the maximum perturbation inthe surface density (defined as the maximum surface densityminus the mean surface density, dashed line). We see thatthe time at which the surface density perturbation comesto exceed the mean surface density, and therefore the pointat which perturbation growth is expected to become non–linear, occurs at ≈ λ with objects of mass m f requires an assumptionabout the radius r of the object in terms of the wavelength.We adopt r = λ/
2, but this likely overestimates the massesof fragments as it includes material whose density is littledifferent from that of the background shell.At later times (but before the formation of the first sinkparticles) the mass function of the p–cores becomes top–heavy (centre and right panels of Figure 7), with a deficit oflow–mass objects, although the slope of the high–mass endof the mass function remains consistent with a power law.In Figure 8 we plot the numbers of p–cores found bythe clumpfinder (before the formation of any sink particles)against the boundedness of the clumps E rat , defined as E rat = | E p | E therm + E k , (7)where E k is the kinetic energy calculated with respect to thecentre of velocity of the p–core, E therm is its thermal energyand E p is the potential energy of the p–core calculated usingthe relative depth of the potential well once the backgroundhas been subtracted. Thus a p-core gravitational potentialdepth is expressed relative to the potential level in its vicin-ity. Values of E rat larger than unity indicate p–cores boundwith respect to their environment, not merely when consid-ered in isolation. Figures 6, 7 and 8 show a clear sequence ofevents in which the evolution of the shell begins to departfrom the predictions of the linear PAGI analysis around thetime when the surface density perturbations come to exceedthe mean shell surface density. The clump mass functionthen ceases to be a power law and the perturbations beginto contract and become gravitationally bound.A fundamental assumption of the thin–shell model isthat all fragments evolve independently of one another in anotherwise unperturbed shell. We used the potential clump–finding code to verify this assumption. The clumpfindingcode is able to follow individual p–cores to see how they gainor lose mass and whether they merge or exchange mass witheach other as time progresses. To make use of this facility,we arbitrarily divide the first 15Myr of the medium–pressureshell’s evolution into ten contiguous 1.5 Myr epochs. In Fig-ure 9 we plot the total number of p–cores extant in eachepoch together with the numbers of clumps formed and de-stroyed in the previous epoch, against time. Initially, thenumber of p–cores destroyed in the previous epoch is al-most equal to the total number of p–cores, implying thatthe objects detected by the clumpfinder are transient enti-ties. However, as the simulation progresses, the total num-ber of p–cores rises and the number destroyed falls, so thatthe objects clearly become more longlived. Around 15Myr,the total number of p–cores levels off (and the numbers ofnew ones formed and of existing ones dissolving declines al-most to zero) implying that the shell evolves to contain astable number of long–lived objects. In Figure 10 we plotthe numbers of p–cores in each 1.5Myr epoch which havesurvived from the previous epoch, the numbers which have c (cid:13) , 1–12 he fragmentation of expanding shells III: Oligarchic accretion and the mass spectrum of fragments . . . N o . o b j ec t s N o . o b j ec t s . . . ⊙ )Mass (M ⊙ )Time: 8.80 Myr . . . N o . o b j ec t s N o . o b j ec t s . . . ⊙ )Mass (M ⊙ )Time: 14.77 Myr . . . N o . o b j ec t s N o . o b j ec t s . . . ⊙ )Mass (M ⊙ )Time: 17.75 Myr . . . N o . o b j ec t s . . . ⊙ )Time: 4.33 Myr . . . N o . o b j ec t s . . . ⊙ ) 7.31 Myr . . . N o . o b j ec t s . . . ⊙ ) 10.29 Myr Figure 7.
Mass functions of clumps identified by the potential–based clump finder in the medium–pressure calculation (top row) andin the high–pressure calculation (bottom row) at three epochs (histograms). Left panels show the mass functions derived from the PAGImodel (solid lines) and centre and right panels show power–laws with a logarithmic slope -2.25 (dashed line). N o . o b j ec t s N o . o b j ec t s .
25 0 . .
75 1 1 .
25 1 . N o . o b j ec t s N o . o b j ec t s .
25 0 . .
75 1 1 .
25 1 . N o . o b j ec t s N o . o b j ec t s .
25 0 . .
75 1 1 .
25 1 . Figure 8.
Boundedness (as defined by Equation 7) of fragments detected by the potential clump–finding algorithm at three epochs inthe medium–pressure calculation. retained more, or less, than half of the particles previouslyassigned to them and the numbers that have merged sincethe previous epoch. The number of p–cores surviving fromthe previous epoch rises rapidly as they become longlived, asopposed to ephemeral. We find that, in the vast majority ofcases an SPH particle assigned to a given p–core remains as-signed to that p–core unless the p–core itself dissolves, andthat mergers or fragmentation of cores are very rare. Ourfindings therefore support the assumption that surface den- sity perturbations do not interact with each other. We alsonote that the total number of p–cores levels off at a value of ∼ ∼
25. Ifthe mass of a fragment is πλ M/ πR , the number of frag-ments with this wavenumber that the shell can support is4 l /π ≈ c (cid:13) , 1–12 James E. Dale, Richard W¨unsch, Rowan J. Smith, Anthony Whitworth, Jan Palouˇs
Time: 22.15 Myr0 . . . N o . o b j ec t s N o . o b j ec t s ⊙ )Mass (M ⊙ ) Time: 23.73 Myr0 . . . N o . o b j ec t s N o . o b j ec t s ⊙ )Mass (M ⊙ ) Time: 26.04 Myr0 . . . N o . o b j ec t s N o . o b j ec t s ⊙ )Mass (M ⊙ ) . . . N o . o b j ec t s ⊙ ) . . . N o . o b j ec t s ⊙ ) . . . N o . o b j ec t s ⊙ ) Figure 11.
Mass functions of sink particles at three epochs in the medium–pressure calculation (top row) and the high–pressurecalculation (bottom row).
Once the p–cores start to become bound, they collapse andform sink particles and the mass function is more easily fol-lowed by tracing the masses of the sinks themselves, shownin Figure 11. We see that the sink–particle mass functionsrapidly come to resemble the clump mass functions pre-sented in the previous section, since there are, by the epochof sink–formation, very few clumps that are not bound andcollapsing. The sink mass functions become ever more top–heavy, and even the high–mass end of the mass function losesits power–law slope. We note that the medium–pressuremass function peaks at a mass of ∼ ⊙ , corresponding tothe mass associated with the most unstable wavenumber of l ∼
25 from Figure 3.Implicit in the thin–shell model is the assumption thata single perturbation will produce a single collapsed objectwhose mass is proportional to the mass of the perturbation,which in turn is given by λ Σ /
16, where λ is the perturba-tion wavelength and Σ is the mean shell surface density atthe time when the fragment collapses. In this picture, theaccretion rate onto each collapsed object thus rises to a peakand falls back down to zero. On the contrary, as shown inFigure 12, we find that, although the accretion rate ontoa given sink is a maximum at the time of the sink parti-cle’s formation and the accretion rate tails off, for most sinkparticles it never falls to zero and they continue accretingfor the duration of the simulation so that, as shown in Fig-ure 13, sink masses are approximately proportional to their ages. We find that for most sink particles, the time–averagedaccretion rate is a strong function of neither age nor mass.The uniformity in the accretion rates of the sink particlesis a consequence of the fact that most of them form in theperiod of time when the shell is almost stationary at its max-imum radius, so that most sinks are born embedded in gasof the same density. Since the shell and the sinks are notin relative motion, the rate of accretion onto a given sink isdetermined only by the local freefall time, which is in turnthe same throughout the shell at any given time.In Figure 14, we show the rate of sink particle formationand the cumulative number of sink particles as functions oftime in the simulation. The number initially rises stronglybefore beginning to tail off around 23 Myr. We showed inFigure 9 that, as the first p–cores become bound, the to-tal number of p–cores stabilises. The final number of sinkparticles in Figure 14 is very similar to the final number ofp–cores, and the tailoff in the rate of sink–formation appearsto be due simply to the fact that there are only a finite num-ber of cores available from which to form sink particles. Thetop–heavy form of the mass function is then a consequenceof the early pulse of sink–particle formation coupled withthe fact that all sink particles share a common environmentand thus accrete at comparable rates. We refer to this pro-cess as oligarchic accretion, since it is the sink particles thatform first that accrete most of the shell’s mass. c (cid:13) , 1–12 he fragmentation of expanding shells III: Oligarchic accretion and the mass spectrum of fragments N o . p - c o r e s N o . p - c o r e s Figure 9.
Total number of p–cores found by the clumpfinderin the medium–pressure run as a function of time (solid line),number of p–cores destroyed in the previous epoch (dashed line)and number of new p–cores formed in the previous epoch (dottedline). Epochs are 1.5Myr in duration..
The purpose of this paper was to simulate the fragmentationof an expanding shell, following the evolution well into thenon–linear regime to see if the predictions of the thin–shelland PAGI models (nominally only valid in the linear regime)can be extrapolated to accurately determine the mass func-tion of fragments formed. We find that the mass functionof objects located by our clump–finding code agrees reason-ably well with the predictions of the thin–shell model in theearly stages of fragmentation before significant numbers offragments become bound. We also find that the assumptionimplicit in the model that fragments do not interact appearsto be sound. However, once p–cores start to become bound,the numbers of p–cores detected levels off at a number con-sistent with the most unstable wavenumber at this epoch.Once the first cores become bound, they suppress the for-mation of new fragments. Once sink particles begin to fromfrom the p–cores, the rate of sink formation is initially high,but tails off as the available cores are consumed. Sink par-ticles continue to accrete at roughly constant rates, sincethey all inhabit a similar environment, so that the initialpulse of sink–particle formation translates into a top–heavymass function in which the mass corresponding to the mostunstable wavelength at the time of the p–cores becomingbound is over–represented with respect to other modes. Theanalysis of Whitworth et al. (1994b), in which the first un-stable wavenumber is taken as representative of the massfunction, describes fragmentation better than the analysisof W¨unsch & Palouˇs (2001) in which all modes contributeto the mass function.We refer to this process as oligarchic accretion, sincethe objects that form first win simply by virtue of beingfirst. This mechanism is very different from the competitiveaccretion described by Bonnell et al. (2001), since in that N o . p - c o r e s N o . p - c o r e s >
50% particlesP-cores with <
50% particlesP-cores merged
Figure 10.
Number of p–cores surviving from previous epoch(solid line), number of p–cores retaining more than half of theparticles assigned to them in the previous epoch (dashed line),number of p–cores retaining less than half of the particles assignedto them in the previous epoch (dotted line) and number of p–coresinvolved in mergers in the previous epoch (dash–dotted line), allplotted as functions of time, and all from the medium–pressurerun. Epochs are 1.5Myr in duration. S i n k m a ss e s ( M ⊙ ) S i n k m a ss e s ( M ⊙ ) . . . . Figure 12.
Evolution of the mass of a randomly chosen 10% ofthe sinks formed in the simulation as a function of time in themedium–pressure calculation. model it is the different environments of the few stars in thedense gas at a cluster centre that enable them to accretemore than their siblings. In the simulations presented here,most of the sink particles form in gas of the same density,around the time when the shell expansion stalls. It is insteadthe time of formation and the non–uniform rate of formationof sink particles which is important in determining the finalmass of a given object. Those objects forming first are able c (cid:13) , 1–12 James E. Dale, Richard W¨unsch, Rowan J. Smith, Anthony Whitworth, Jan Palouˇs S i n k m a ss e s ( M ⊙ ) S i n k m a ss e s ( M ⊙ ) . . . Figure 13.
Masses of sink particles as a function of their age asmeasured at a time of 31 Myr in the medium–pressure calculation. S i n k f o r m a t i o n r a t e ( a r b i t r a r y un i t s ) S i n k f o r m a t i o n r a t e ( a r b i t r a r y un i t s ) . . . N o . s i n k s N o . s i n k s Figure 14.
Rate of formation (solid line) and total number(dashed line) of sink particles in the medium–pressure run plottedagainst time . to accrete more mass. The mass function becomes skewedbecause the rate of sink particle formation declines as thesupply of p–cores formed during the linear phase of the shellsevolution is consumed. Klessen & Burkert (2000) observe asimilar process in the evolution of a turbulent protocluster,where the mass function becomes skewed towards highermasses as the star formation efficiency becomes large andthe gas reservoirs required to form new cores are depleted.The accretion process may also be aided by the fact that, ata time of ≈
20 Myr, the shell begins to contract. However,the contraction is slow compared to the rate of sink particleformation and accretion, and the fragment mass function isalready strongly top–heavy at this epoch, so contraction ofthe shell cannot be the primary driver of oligarchic accre-tion. It is possible that, in the case of a shell sweeping up massas it expands, so that reserves of fresh gas in the shell areconstantly replenished, p–cores would be able to form at alltimes, resulting in a mass function more closely resemblinga power law. However, such a shell would experience theVishniac instability (Vishniac 1983), probably altering themass spectrum of fragments in the non–linear portion of theshell’s evolution. This is important, since our results suggestthat the most unstable mode at the time when bound ob-jects begin to form is over–represented in the mass function.To unequivocally demonstrate that accretion of fresh mate-rial would allow unabated p–core formation and produce apower–law mass function would require that the shell’s gasreservoir be replenished in some way during sink formation.This is difficult in the case presented here, since most of thesinks form when the shell is almost stationary. The only wayto replenish the gas would then be to actively feed matterinto the shell, perhaps allowing it to infall from a reservoirjust outside the shell’s maximum radius, a highly artificialconstruction. We therefore defer answering this question toa later paper in which we analyse the effect of the Vishniacinstability on a momentum–driven shell sweeping up an ex-ternal medium.The mass functions produced by our simulations bearlittle resemblance to any known stellar or cluster mass func-tion. As explained in Section 3, the mass resolution of oursink particles is good enough to resolve low–mass stars, butthe linear resolution implied by the accretion radii (0.05pc)is rather large. Our sink particles should strictly be regardedas binaries or small multiple systems so that our mass func-tions, in the terminology of ? , are Multiple Star Mass Func-tions (MSMFs), not Single Star Mass Functions (SSMFs).This would be true to an even greater degree if our resultswere rescaled to higher masses. ? find, not surprisingly, thatthe SSMF contains fewer high-mass objects and more low–mass objects relative to the MSMF, and that the SSMFpeaks at lower masses. The SSMF from our simulationswould probably therefore be slightly less deficient in low–mass objects and contain fewer high–mass objects than themass functions we present, but the effects of multiplicity willnot quantitatively change our finding that the mass func-tions become very top–heavy due to oligarchic accretion.Dawson et al. (2008) observed the mass function ofCO clumps in the Carina Flare supershell and observedthat their mass spectrum follows a power–law. Althoughthey caution against over–interpretation of their data, theirclump mass measurements imply that the clumps are notgravitationally bound (in the sense that their inferredmasses are below their virial masses). This result is con-sistent with our findings that the as–yet–unbound objectsproduced in the linear phase of shell fragmentation havea power law mass function (although note that the powerlaw we derive is considerably steeper than that inferred byDawson et al. (2008)). Our findings suggest, however, thatthe mass function of these unbound fragments cannot nec-essarily be used to predict the mass function of stars or starclusters that will eventually form. There are as yet no obser-vations of the stellar or cluster mass functions generated bythe fragmentation of expanding shells, owing to the extremedifficulty of making such a measurement.The initial conditions used in our simulations, that ofa shell which is pressure–confined and yet does not sweep c (cid:13) , 1–12 he fragmentation of expanding shells III: Oligarchic accretion and the mass spectrum of fragments up any mass are somewhat unusual and the peculiar massfunction we obtain implies that most stars or clusters donot form in shells of this type. However, such shells mustsurely occur sometimes, since a shell may expand into thehot rarefied ISM and be pressure–confined while accretingnegligible mass. Our results demonstrate that shells of thisnature will produce anomalously high numbers of massivestars or clusters and may therefore lead to the propagationof star formation by triggering. We find that, in the case of a momentum–driven shellpressure–confined internally and externally, the PAGI modelfaithfully reproduces the mass function of fragments pro-duced by the gravitational instability during the time forwhich perturbations in the shell surface density are rela-tively small and the behaviour is linear.In the case studied here of a shell which does not sweepup mass, we find that a form of accretion which we termoligarchic accretion operates within the shell whereby thefirst objects to form become the most massive because theyaccrete from the gas reservoir for longer. Coupled with thefinite number of p–cores formed during the linear fragmenta-tion of the shell, this leads to a very top–heavy mass functionin which the most unstable wavelength at the time when per-turbations first become bound is over–represented. A shellwhose gas reservoir is replenished may form fragments con-tinuously, resulting in a mass function more closely resem-bling a power law.Since the mass function generated by our simulationsis so unusual, we infer that shells of the type we have stud-ied cannot make a strong contribution to the global stellarpopulation, although no direct measurements of the massfunctions produced by expanding shells exist against whichwe could test our conclusions.
The authors thank the referee, Sven Van Loo, for inci-sive and helpful comments which considerably improved thestructure of the paper. JED acknowledges support from aMarie Curie fellowship as part of the European Commis-sion FP6 Research Training Network ‘Constellation’ undercontract MRTN–CT–2006–035890. JED, RW and JP ac-knowledge support from the Institutional Research PlanAV0Z10030501 of the Academy of Sciences of the CzechRepublic and project LC06014–Centre for Theoretical As-trophysics of the Ministry of Education, Youth and Sportsof the Czech Republic. AW gratefully acknowledges the sup-port of grant ST/H001530/1 from the UK Science and Tech-nology Facilities Council.
The thin shell dispersion relation ω ( l ) = − AVR + r BV R − c s l R + 2 πG Σ lR (8) may be recast in a dimensionless form with the help of anexpansion law – a prescription for the variations in the shellradius R , expansion velocity V and surface density Σ withtime. The shells described in this work are ballistic and offixed mass M , expanding due to their inertia into a vacuumand decelerated by self–gravity alone.We define a dimensionless radius ξ as ξ = RR max = 1 − K | W | = 1 − V RGM (9)where K = MV / W = ( GM ) / (2 R ) are the kineticand potential energy respectively, and R max = R − K / | W | (10)is the maximum radius to which the shell expands (the initialtotal energy K + W at radius R is assumed to be negativeso that the shell is bound). The free fall time obtained fromKepler’s third law is t ff = πR / GM ) / (11)and the most unstable dimensionless wavenumber at R max is l max = πG Σ R max c s = GM c s R max . (12)Inserting (9) and (12) into (8) and normalizing it by (11) weobtain w ( l, ξ, l max ) = ω ( l, R, V, Σ , c s ) t ff = − Aπ − ξ ) / ξ − / + π × (cid:20) B (1 − ξ ) ξ − − l l max ξ − + 12 lξ − (cid:21) / . (13)Using numerical factors A = 3 / B = 1 / w ( l, ξ, l max ) = − π − ξ ) / ξ − / + π (cid:20) ξ − (1 + 2 l ) − ξ − (cid:18) l l max (cid:19)(cid:21) / . (14)The degree of fragmentation at a given wavenumber l isdescribed by the fragmentation integral I f ( l, l max ) = Z tt ω ( l, t ′ , l max ) dt ′ = Z ξξ w ( l, ξ ′ , l max ) t ff dtdξ dξ ′ . (15)Since dξdt = 1 R max dRdt = VR max (16)and VR max = (1 − ξ ) / (cid:18) GMR (cid:19) / R max = (1 − ξ ) / (cid:18) GMR max (cid:19) / R max (cid:18) R max R (cid:19) / = (cid:18) − ξξ (cid:19) / π t ff , (17) c (cid:13) , 1–12 James E. Dale, Richard W¨unsch, Rowan J. Smith, Anthony Whitworth, Jan Palouˇs where equations (9) and (11) have been used, we may write dtdξ = (cid:18) ξ − ξ (cid:19) / t ff π . (18)Inserting (18) into (15) yields the fragmentation integral ina form I f ( l, l max ) = 2 π Z ξξ w ( l, ξ ′ , l max ) (cid:18) ξ ′ − ξ ′ (cid:19) / dξ ′ = − Z ξξ ξ ′− dξ ′ + 12 × Z ξξ (cid:20) lξ ′ (1 − ξ ′ ) − l /l max ξ ′ (1 − ξ ′ ) (cid:21) / dξ ′ . (19)Using the Sage symbolic manipulation package (available at ), the fragmentation integral becomes I f ( l, l max ) = −
32 ln (cid:18) ξξ (cid:19) + 12 (2 l + 1) / × ln (cid:26) ξξ l + 1) l max ] / α ( ξ ) + βξ + (4 l + 2) l max l + 1) l max ] / α ( ξ ) + βξ + (4 l + 2) l max (cid:27) +12 (cid:18) l max + l l max (cid:19) / × ln (cid:26) l max + l ) / α ( ξ ) + (2 l max + 2 l ) ξ + β l max + l ) / α ( ξ ) + (2 l max + 2 l ) ξ + β (cid:27) , (20)where α ( x ) = [( l max + l ) x + βx + (2 l + 1) l max ] / (21)and β = ( − l − l max − l . (22) REFERENCES
Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277,362Benz W., 1990, in Buchler J. R., ed., Numerical Modellingof Nonlinear Stellar Pulsations Problems and ProspectsSmooth Particle Hydrodynamics - a Review. pp 269–+Bonnell I. A., Bate M. R., Clarke C. J., Pringle J. E., 2001,MNRAS, 323, 785Churchwell E., Povich M. S., Allen D., Taylor M. G., MeadeM. R., Babler B. L., Indebetouw R., Watson C., Whit-ney B. A., Wolfire M. G., Bania T. M., Benjamin R. A.,Clemens D. P., Cohen M., Cyganowski C. J., JacksonJ. M., Kobulnicky H. A., Mathis J. S., 2006, ApJ, 649,759Churchwell E., Watson D. F., Povich M. S., Taylor M. G.,Babler B. L., Meade M. R., Benjamin R. A., IndebetouwR., Whitney B. A., 2007, ApJ, 670, 428Dale J. E., Bonnell I. A., Whitworth A. P., 2007, MNRAS,375, 1291Dale J. E., W¨unsch R., Whitworth A., Palouˇs J., 2009,MNRAS, 398, 1537Dawson J. R., Kawamura A., Mizuno N., Onishi T., FukuiY., 2008, PASJ, 60, 1297Deharveng L., Lefloch B., Kurtz S., Nadeau D., Pomar`esM., Caplan J., Zavagno A., 2008, A&A, 482, 585Deharveng L., Lefloch B., Massi F., Brand J., Kurtz S.,Zavagno A., Caplan J., 2006, A&A, 458, 191 Ehlerov´a S., Palouˇs J., 2005, A&A, 437, 101Elmegreen B. G., 1994, ApJ, 427, 384Elmegreen B. G., Lada C. J., 1977, ApJ, 214, 725Klessen R. S., Burkert A., 2000, ApJS, 128, 287K¨onyves V., Kiss C., Mo´or A., Kiss Z. T., T´oth L. V., 2007,A&A, 463, 1227Monaghan J. J., 1992, ARA&A, 30, 543Rod´on J. A., Zavagno A., Baluteau J., Anderson L. D.,Polehampton E., Abergel A., Motte F., Bontemps S., AdeP., Andr´e P., 2010, ArXiv e-printsSmith R. J., Clark P. C., Bonnell I. A., 2009, MNRAS, 396,830Vishniac E. T., 1983, ApJ, 274, 152Whitworth A. P., Bhattal A. S., Chapman S. J., DisneyM. J., Turner J. A., 1994a, A&A, 290, 421Whitworth A. P., Bhattal A. S., Chapman S. J., DisneyM. J., Turner J. A., 1994b, MNRAS, 268, 291Williams J. P., de Geus E. J., Blitz L., 1994, ApJ, 428, 693Wunsch R., Dale J. E., Palous J., Whitworth A. P., 2010,ArXiv e-printsW¨unsch R., Palouˇs J., 2001, A&A, 374, 746Zavagno A., Deharveng L., Comer´on F., Brand J., MassiF., Caplan J., Russeil D., 2006, A&A, 446, 171Zavagno A., Russeil D., Motte F., Anderson L. D., De-harveng L., Rodon J. A., Bontemps S., Abergel A., Ba-luteau J., Sauvage M., Andr´e P., Hill T., White G. J.,2010, ArXiv e-prints c (cid:13)000