Fragmentation of vertically stratified gaseous layers: monolithic or coalescence-driven collapse

8.5 MB·

by  František Dinnbier

Mon.Not.R.Astron.Soc.000,1–21(2017) Printed12January2017 (MNLATEXstylefilev2.2)
Fragmentation of vertically stratified gaseous layers:
monolithic or coalescence–driven collapse
Frantiˇsek Dinnbier1,2(cid:63), Richard Wu¨nsch1, Anthony P. Whitworth3 and Jan Palouˇs1
1Astronomical Institute, Academy of Sciences of the Czech Republic, Boˇcn´ı II 1401, 141 00 Prague, Czech Republic
7 2Charles University in Prague, Faculty of Mathematics and Physics, V Holeˇsoviˇck´ach 2, 180 00 Prague, Czech Republic
1 3School of Physics and Astronomy, Cardiff University, Cardiff CF24 3AA, Wales, UK
0
2

n Accepted2016December21.Received2016December20;inoriginalform2016August09
a
J

0 ABSTRACT
1 We investigate, using 3D hydrodynamic simulations, the fragmentation of pressure-
confined, vertically stratified, self-gravitating gaseous layers. The confining pressure

] is either thermal pressure acting on both surfaces, or thermal pressure acting on one
A
surface and ram-pressure on the other. In the linear regime of fragmentation, the dis-
G persion relation we obtain agrees well with that derived by Elmegreen & Elmegreen
. (1978),andconsequentlydeviatesfromthedispersionrelationsbasedonthethinshell
h
approximation (Vishniac 1983) or pressure assisted gravitational instability (Wu¨nsch
p
et al. 2010). In the non-linear regime, the relative importance of the confining pres-
-
o sure to the self-gravity is a crucial parameter controlling the qualitative course of
r fragmentation.Whenconfinementofthelayerisdominatedbyexternalpressure,self-
t
s gravitatingcondensationsaredeliveredbyatwo-stageprocess:firstthelayerfragments
a
into gravitationally bound but stable clumps, and then these clumps coalesce until
[
they assemble enough mass to collapse. In contrast, when external pressure makes a

1 smallcontributiontoconfinementofthelayer,thelayerfragmentsmonolithicallyinto
v gravitationally unstable clumps and there is no coalescence. This dichotomy persists
8 whether the external pressure is thermal or ram. We apply these results to fragments
7 forminginashellsweptupbyanexpandingHiiregion,andfindthat,unlesstheswept
7
up gas is quite hot or the surrounding medium has low density, the fragments have
2 low-mass(<3M ),andthereforetheyareunlikelytospawnstarsthataresufficiently
0 ∼ (cid:12)
massive to promote sequential self-propagating star formation.
.
1
0 Key words: stars: formation – ISM: H ii regions – ISM: kinematics and dynam-
7 ics – Physical processes: instabilities – Physical processes: hydrodynamics – Physical
1 processes: waves
:
v
i
X
r
a 1 INTRODUCTION propagates itself sequentially. Otherwise, star formation is
quenched. Various triggering mechanisms are discussed in
Massive stars (M >8M ) strongly influence the interstel-
(cid:63)∼ (cid:12) Elmegreen (1998).
lar medium (ISM) surrounding them, mainly via photoion-
isation, stellar winds and supernova explosions. Elmegreen From an observational perspective, shells are common
&Lada(1977)proposeamechanism(Collect and Collapse) morphological structures in the ISM. Ehlerova´ & Palouˇs
wherebyanover-pressuredHiiregion,drivenbyyoungmas- (2005)havedetectedmorethan600shellsinHi.Churchwell
sivestars,expandsintodensemoleculargas.Theexpansion etal.(2007)haveidentified322completeorpartialringsin
induces a spherical shock, and the surrounding gas accu- the infrared. Simpson et al. (2012) list more than 5000 in-
mulates between the shock and the ionisation front. The fraredshells.Deharvengetal.(2010)findthatatleast86%
resulting shell of cool gas increases in mass, and eventually of the infrared shells identified by Churchwell et al. (2007)
fragmentstoformanewgenerationofstars.Acrucialissue encircle H ii regions ionised by O- and early B-type stars,
is the maximum mass of these newly formed stars. If some suggesting that the shells are due to feedback from these
of them are sufficiently massive to excite a new H ii re- stars. After a careful examination, Deharveng et al. (2005)
gion, the process can repeat recursively and star formation find17infraredshellsthatarecandidatesfortheCollectand
Collapse mechanism. Evidence for propagating star forma-
tion has also been reported in significantly larger H i shells
(cid:63) E-mail:[email protected] (Dawson et al. 2008, 2011; Egorov et al. 2014), indicating
(cid:13)c 2017RAS
2 F. Dinnbier et al.
that feedback operates over a large range of scales. Further gaseous layers, in both the low and high ambient pressure
observationalsupportforpropagatingstarformationcomes cases,andtocomparetheresultswiththeanalyticorsemi-
from age sequences of OB associations and star clusters in analyticestimatesderivedinpreviousstudies.Thiscompar-
thevicinityofstarformingregionse.g.(Blaauw1964,1991; ison addresses three issues: (i) dispersion relations in the
Brown et al. 1994; Bik et al. 2010) linear regime of fragmentation; (ii) the elapsed time before
In order to estimate the properties of fragments con- thelayerformsgravitationallyboundfragments,andthere-
densing out of a swept up shell, Elmegreen & Elmegreen sulting fragment masses; and (iii) the possibility that mode
(1978, hereafter E78), Doroshkevich (1980) and Lubow & interaction leads to fragments distributed on a regular pe-
Pringle (1993) investigate the stability of a vertically re- riodically repeating pattern. We study only a small square
solvedself–gravitatinglayer1 ,andderiveasemi-analytical patchonthelayer,sothatwecanattaingoodresolutionin
formula for the corresponding dispersion relation. Vishniac directionsperpendiculartothelayer,evenwhenthelayeris
(1983, hereafter V83) derives the dispersion relation for a significantlycompressed.Inadditiontolayersconfinedfrom
self–gravitating infinitesimally thin shell, and this has been both sides by thermal pressure, we also simulate layers ac-
usedtoestimatethepropertiesofclumpscondensingoutof cretingontoonesurface.Weusetheseresultstoinvestigate
fragmentingshells.Whitworthetal.(1994a)andElmegreen the fragmentation of a shell driven by an expanding H ii
(1994) analyse the fragmentation of accreting shells, while region — using an analytic solution to account for the ex-
Wu¨nsch & Palouˇs (2001) consider non-accreting shells ex- pansion of the H ii region — and contrast our results with
panding into a vacuum. thoseobtainedanalyticallybyWhitworthetal.(1994a),and
Thesedispersionrelationsarebasedonlinearperturba- by Iwasaki et al. (2011b) who performed simulations with
tiontheory,sotheyarerelevantonlyaslongastheperturb- small perturbing amplitudes.
ing amplitudes are small. Miyama et al. (1987a,b) extend Thepaperisorganisedasfollows.Section2reviewsthe
the linearised solution of Elmegreen & Elmegreen (1978) mostimportantpropertiesofself–gravitatinglayersandthe
intothenon-linearregime,byincludingsecondorderterms. analyticalestimates(E78,V83,W10)ofthedispersionrela-
Theyfindthatalayerbreaksintofilamentsthatbecomein- tions.Theseestimatesareusedforcomparisonwiththenu-
creasingly slender with time. In contrast, Fuchs (1996) pro- mericalresults.Section3describesappliednumericalmeth-
poses that, in the non-linear regime, the fragments form a odsandinitialconditions.Section4describessimulationsof
semi-regular hexagonal pattern on the surface of the layer. layers confined on both sides by thermal pressure (due to a
In order to test the validity of these analytic deriva- veryhotrarefiedgas),andseededwithmonochromaticper-
tions, Dale et al. (2009) obtain a dispersion relation, based turbations.Section5describessimulationsoflayersconfined
onnumericalsimulationsofexpanding,non-accretingshells, on both sides by thermal pressure, and seeded with poly-
confinedbyconstantexternalpressure.Theyfindsignificant chromatic perturbations. Section 6 describes simulations in
differencesbetweentheirdispersionrelationandthatofV83 which we explore whether fragments forming in layers tend
(basedonthethinshellapproximation).Insubsequentwork, tobearrangedintoregularpatterns.Section7describessim-
Wu¨nsch et al. (2010, hereafter W10) explain the difference ulations of layers in which one side is confined by thermal
by modifying V83 to include the effect of pressure confine- pressure and the other by the ram pressure of a homoge-
ment. They call the mechanism underlying their dispersion neous plane-parallel inflow. Section 8 considers an H ii re-
relation Pressure Assisted Gravitational Instability (PAGI gionexpandingintoahomogeneousmedium,andestimates
Wu¨nsch et al. 2010). Since Dale et al. (2009) simulate frag- the time at which the swept up shell fragments, and the
mentation of a whole shell, and the shell gets thinner with properties of the fragments. Section 9 discusses the results,
increasing external pressure, they are only able to resolve andSection10summarisesourmainconclusions.Appendix
shells which are confined by low to moderate values of the Adescribesouralgorithmdevelopedforfindinggravitation-
externalpressure,andthereforetheW10dispersionrelation ally bound fragments.
is unverified in the case of high external pressure. However,
it is in the limit of high external pressure that W10 dif-
fers substantially from E78. Van Loo et al. (2014) investi-
gateself-gravitatinglayerspermeatedbymagneticfieldsand 2 LINEARISED THEORY OF LAYER
present a non-magnetic control run. Although their results FRAGMENTATION
differ from V83, they are in agreement with both E78 and
2.1 The unperturbed state
W10, because their layer is confined by a very low external
pressure. Iwasaki et al. (2011b) also derive a numerical dis- Before introducing the dispersion relations describing the
persion relation that differs significantly from that of V83. growthrateofperturbationsintheself-gravitatingpresure-
Our main aims are to simulate the fragmentation of confined layer, we review the unperturbed configuration.
The model is also used to generate the initial conditions
for our simulations.
1 We distinguish a layer, which ideally is plane-parallel strati- We assume that the layer is initially in plane-parallel
fied,fromashell,whichideallyispart,orthewhole,ofahollow stratified hydrostatic equilibrium, i.e. it is infinite in the x
spherically symmetric structure. Although most of the observa- and y directions, its normal points in the z direction, and
tions motivating this work concern shells, the configurations we allquantitiesarefunctionsofz only.Thelayerisisothermal
simulate in this paper are layers. Then we apply results derived
withsoundspeedc ,sothedensitydistributionisasderived
for layers in the analysis of star formation triggered by expand- S
bySpitzer(1942)andGoldreich&Lynden-Bell(1965),i.e.
ingHiiregions—sointhatanalysiswealsousethetermshell,
becausethespatialcurvatureandvelocitydivergenceoftheshell ρ
ρ(z)= O , (1)
canbeneglected. cosh2(z/H )
O
(cid:13)c 2017RAS,MNRAS000,1–21
Fragmentation of stratified gaseous layers 3
where ρ is the density on the midplane (z=0), rateω isattainedforwavenumberk whichisalways
O FAST FAST
HO = (cid:112)2πcSGρ , (2) kapMpArXoxainmdakteFlAySTkMaAreX/d2e.noTthede λwMavAeXleanngdthλsFAcoSTrr,ersepsopnedcitnivgeltyo.
O Whendescribingaparticularanalyticalestimateofthedis-
is the vertical scale height, and G is the gravitational con- persion relation, the subscript begins with its name, e.g.
stant.Thesurfacesofthelayerareatz=±z ,andoutside k is the wavelength k for E78. The abbrevia-
MAX E78.FAST FAST
this (|z|>z ) there is a hot gas with negligible density tions E78, V83 and W10 refer to the particular dispersion
MAX
which exerts an external pressure P . relation, not to the paper where they are firstly described.
EXT
The layer is fully characterised by P , c and its The unstable branches of the dispersion relations for both
EXT S
surface density Σ . A dimensionless parameter constructed self–gravity (A = 0.99) and external pressure dominated
O
from these quantities (Elmegreen & Elmegreen 1978), (A=0.18) layers are shown in the right column of Fig. 2.
1
A= , (3)
(cid:113)
1+(cid:0)2P /πGΣ2(cid:1) 2.2.1 Dispersion relation for the thin shell approximation
EXT O
To estimate the dispersion relation of an expanding shell,
reflects the relative importance of self-gravity and external
Vishniac (1983) reduce the problem to two dimensions by
pressureinholdingthelayertogether.ParameterAincreases
integrating the continuity, Euler, and Poisson’s equation
monotonically from nearly 0 (external pressure dominated;
P (cid:29)GΣ2) to 1 (self–gravity dominated). through the thickness of the shell. In planar geometry, the
EXT O thin shell dispersion relation becomes
From pressure equilibrium at the surfaces, P =
EXT
ρ(z )c2, it follows ω2 (k)=2πGΣ k−c2k2. (7)
MAX S V83 O S
(cid:18)ρ c2(cid:19) (cid:18) 1 (cid:19) The stability of modes is determined by the imbalance be-
z =H arcosh−1 O S =H arcosh−1 √ . tween self-gravity and internal pressure gradient, and there
MAX O P O 1−A2
EXT is no contribution from external pressure.
(4)
Half–thickness of a stratified layer is defined by H ≡
HT 2.2.2 Dispersion relation for pressure assisted
Σ /(2ρ ). Equations (3), (1) and (2) then yield
O O gravitational instability
Σ c2 c2A2
H = O S = S , (5) In order to evaluate the influence of external pressure,
HT 2P +πGΣ2 πGΣ Wu¨nschetal.(2010)investigateperturbationswiththeform
EXT O O
of an oblate spheroid, embedded in a layer. The spheroid is
and the midplane density, ρ , depends on Σ as
O O homogeneous and confined by external pressure P . Its
2P +πGΣ2 πGΣ2 semi-major axis is r, its semi-minor axis is the layeErX’sThalf-
ρ = EXT O = O. (6)
O 2c2 2c2A2 thickness, H , and its total mass is M. Radial excursion
S S HT
ofanelementontheequatorofthespheroidisregulatedby
Pressure dominated layers are of almost uniform density
the equation of motion
(z (cid:39) H ), while self–gravity dominated layers have
a MprAoXnounceHdTdensity maximum on the midplane (z (cid:29) 3GM (cid:26) cos−1(H /r) H /r (cid:27)
MAX r¨ = − HT − HT
H ). 2r2 (1−(H /r)2)3/2 1−(H /r)2
HT HT HT
20πrH P 5c2
− HT EXT + S (8)
3M r
2.2 Analytical estimates of the dispersion relation
(Boyd & Whitworth 2005).
Inthissection,wereviewandcomparetheassumptionsun-
Wu¨nschetal.(2010)equatetheinstabilitygrowthrate
derlying the dispersion relations derived by E78, V83 and
totheradialcontractionrateofthespheroid.Collapsefrom
W10. To simplify the discussion, we neglect the effects of
radius r , by a small factor (cid:15) to radius (1−(cid:15))r , in time
shell curvature and velocity divergence, by setting the shell t , fulfilsO((cid:15)−1)r = r˙ t + 1r¨ t2. Using ω =O1/t and
radius to infinity. The results are therefore applicable to a (cid:15) O O (cid:15) 2 O (cid:15) (cid:15)
r=λ/2=π/k, this yields
plane–parallel layer, and can also be applied to a shell, as
longastheunstablewavelengthsaremuchshorterthanthe ω2 (k) = − r¨O
radius of the shell. W10 2(cid:15)r
O
The dispersion relation gives the perturbation growth 1(cid:26) 5c2k2
= − S
rate ω as a function of wavenumber k. In planar geome- (cid:15) 2π2
tpruyr,ealyndrefaolrotrhpeudrieslpyeirmsiaogninraelrayt.ioInnsthinisqpuaepsetrio,nw,eωadisopetiththeer +3GΣOk(cid:32)cos−1(kHHT/π)−(cid:112)k2HH2T/π2−1(cid:33)
convention that growing instability is described by the pos- 4 (1−k2H2 /π2)3/2
HT
itive real part of ω, which we denote for simplicity ω, and 10P c2k2 (cid:27)
we omit the oscillating imaginary part. The characteristic + EXT S . (9)
3π2(2P +πGΣ2)
timescale of perturbation growth, its e–folding time, is de- EXT O
fined as t =1/ω. The terms on the right hand side of Eq. (9) represent the
EFOLD
Forallthedispersionrelationsinquestion,thereisafi- gradientofinternalpressure(secondline),self-gravity(third
nite range of unstable wavenumbers (0,k ), where k line), and confinement by external pressure (fourth line).
MAX MAX
is the highest unstable wavenumber. The maximum growth Wu¨nsch et al. (2010) suggest setting (cid:15)∼0.1, but, although
(cid:13)c 2017RAS,MNRAS000,1–21
4 F. Dinnbier et al.
Figure 1.ComparisonbetweenthedispersionrelationsderivedbyE78,V83andW10,asafunctionofA. Leftpanel:themarginally
stablewavenumberk ,forE78andW10,normalisedtok .Rightpanel:Massofthefragmentwiththehighestgrowthrate,
MAX V83.MAX
M ,normalisedtothemidplaneJeansmass,M .Fragmentsformedwithmassbelow∼M (thinline)aregravitationallystable.
FAST J J
(cid:15) affects the magnitude of ω , it does not influence the 2.2.4 Comparison between the analytical estimates of the
W10
range of unstable wave-numbers. dispersion relation
2.2.3 Dispersion relation for a vertically stratified layer TheleftpanelofFigure1showsthedependenceofk nor-
MAX
malisedtok onparameterAfordifferentdispersion
Elmegreen & Elmegreen (1978) obtain the dispersion rela- V83.MAX
relations. For self–gravity dominated layers, both W10 and
tion for a self-gravitating, pressure confined, semi-inifinite
E78 predict almost identical ranges of unstable wavenum-
layerinhydrostaticequilibriumbysolvingthesystemofper-
bers (k /k → 0.966 as A → 1), but V83
turbed continuity, Euler and Poisson equations. Kim et al. E78.MAX W10.MAX
predicts a broader range with k /k →0.518.
(2012)revisittheirstudyandofferanadditionalinsight.Ac- W10.MAX V83.MAX
With increasing external pressure (decreasing A), the dif-
cording to Kim et al. (2012), the dispersion relation can be
ference between E78 and V83 diminishes. However, when
written A<0.5,E78extendstomuchhigherwavenumbers,indicat-
∼
ing that shorter wavelengths are unstable. And, as A → 0,
ω2 =2πGΣ k(ηF +(1−η)F )−c2 k2. (10) W10 and E78 have different limiting behaviour: k
E78 O J D EFF W10.MAX
tends to a constant (k /k →2.216), whereas
W10.MAX V83.MAX
Here c is the effective sound speed. F and F are re- k ∝H−1 ∝A−2.
EFF J D E78.MAX HT
duction factors for self-gravity. Parameter η is the fraction
oftheperturbedsurfacedensitythatisattributabletocom-
pressionofmaterialnearthecentreofaproto-fragment(and
hence near the mid-plane of the layer). The rest of the per-
Sincethewavelengthλ hasthehighestgrowthrate,
turbedsurfacedensityisduetocorrugationsonthesurface FAST
we assumethatthe fragmentmass M isapproximately
ofthelayer.Thusηcontrolstherelativeimportanceofcom- FAST
the mass confined inside a circle of radius λ /2, i.e.
pressionalandsurface-gravitywaves.Theformerareimpor- FAST
M =πΣ (λ /2)2 =πΣ (π/k )2.MassesM
tant in self–gravity dominated layers, while the latter are FAST O FAST O FAST FAST
normalised to the midplane Jeans mass M are plotted in
importantinpressuredominatedlayers.c ,η,F andF J
EFF J D the right panel of Fig. 1. V83 and W10 predict that when
are complicated functions of k, which can not generally be
formed, the fragments are already gravitationally unstable
expressed in closed form.
(M /M >1) for any A. On the other hand, E78 pre-
The growth rate ωE78.FAST of the most unstable dictFsASgTravitJa∼tionally unstable fragments only when A>0.5.
wavenumberk inthelimitA→0(layersdominated ∼
E78.FAST According to E78, the layer breaks into gravitationally sta-
by external pressure) is (eq. (47) in Kim et al. (2012))
ble fragments for lower values of A, thus predicting quali-
tatively different scenario than V83 and W10. We simulate
(cid:18)πGΣ (cid:19)2
ω2 =0.276×2πGρ =0.276 O . (11) and discuss evolution of pressure dominated layers (low A)
E78.FAST O c A in Sections 5.3 and 9.1.
S
(cid:13)c 2017RAS,MNRAS000,1–21
Fragmentation of stratified gaseous layers 5
3 METHODS AND INITIAL CONDITIONS Klessen 1997). By computing the appropriate limit of the
standard Ewald method, we find formulae for the gravita-
3.1 Numerics
tional acceleration and potential in closed form for a con-
Allthehydrodynamicsimulationspresentedinthepaperare figuration with mixed BCs; the derivation is described in
performedwiththeMPI-parallelisedflash4.0code(Fryxell (Wu¨nsch et al. in preparation). The modification is used to
et al. 2000). flash is an amr code based on the paramesh calculatethegravitationalfieldinthesimulationspresented
library(MacNeiceetal.2000).Thehydrodynamicequations here.
are solved by the piecewise parabolic method (Colella & The hydrodynamic BCs are periodic in the x and y
Woodward 1984). directions, and reflecting in the z direction for all but ac-
Self-gravity is calculated using an octal tree code creting simulations. For accreting simulations, the BCs are
(Wu¨nschetal.inpreparation)whichoffersthreeacceptance inflowfromthetopofthecomputationaldomain,anddiode
criteria for interaction between the target point (where the from the bottom to prevent reflections of waves.
accelerationisevaluated)andanode(thesourceofthegrav- The column density needed for cooling the warm am-
itational force): bient gas intermixed with the cold layer during accreting
runs (see Section 3.3 for details) is calculated using module
(i) The algorithm invented by Barnes & Hut (1986),
TreeRay/OpticalDepth of the tree code.
which accepts nodes seen from the target point at an an-
Thegridcellsarecubicinallthesimulations.Thehalf-
gle smaller than a specified value. This algorithm is purely
thickness of the layer (H ; see Eqn. 5 below), is always
geometricsinceitdoesnottakeintoaccountdistributionof HT
at least four grid cells (see Tables 1 to 5). Consequently,
mass inside nodes, or their relative contribution to the net
the Jeans length is always resolved by more than four grid
gravitational force.
cells, and the simulations satisfy the criterion for avoiding
(ii) A set of criteria using the node size, mass and op-
artificial gravitational fragmentation, as given by Truelove
tionally higher multipole moments of the mass distribution
et al. (1997). We have performed successful convergence
withinthenode.Theyeitherestimatetheupperlimitonthe
teststhroughouttherangeofconditionssimulated.Wehave
acceleration error of the cell-node interaction ∆amax from
alsocheckedthattheportionofthelayerinsidethecompu-
Equation 9 of Salmon & Warren (1994), or they use the
tational domain is sufficiently large, i.e. that the periodic
approximation by Springel (2005),
copies in the x and y directions do not significantly influ-
GM (cid:18)h(cid:19)p+1 ence the properties of fragments.
∆amax = , (12)
(p) d2 d
where M and h are the node mass and size, respectively, d 3.2 Initial conditions for the layer
is the distance between the cell and the node mass centre,
Eachmodelstartswithalayerwithpropertiesdescribedin
and p is the order of the multipole expansion.
Section2.1seededwithaperturbation.Theperturbationis
(iii) Anexperimentalimplementationofthe”sumsquare”
eitherasinglemodegivenbytheeigenfunctionforacoustic-
criterion of Salmon & Warren (1994), which controls the
surface-gravitymodesorwhitenoise.Theexactformofthe
sumofallerrorsoriginationfromallindividualnodecontri-
perturbation is described at the beginning of a correspond-
butions.
ing Section (4 to 7).
We invoke the criterion (ii) with Equation (12) and p = 1 Particular values for the layer parameters are adopted
in all our simulations because our tests show that it pro- fromIwasakietal.(2011b)whoinvestigateanHiiregionex-
vides the best compromise between the code performance citedbya41M star,expandingintoamediumwithnum-
(cid:12)
and accuracy. berdensityn=103cm−3.Attimet=0.81Myr,theshellhas
We convert dense gaseous condensations into sink par- radiusR=3.86pc,andsurfacedensityΣ =0.0068gcm−2.
O
ticles according to conditions described in Federrath et al. We use this value of Σ for the initial conditions of all the
O
(2010). To create a sink particle inside a particular cell, all simulations presented here, and vary A by changing P .
EXT
following conditions must be fulfilled: However, we note that the model is sufficiently simple that
its physical parameters, (Σ ,c ,P ), can be rescaled ar-
(i) the cell is at the highest refinement level, O S EXT
bitrarily, as long as A is unchanged, i.e. P ∝Σ2.
(ii) the cell contains minimum of gravitational potential, EXT O
The sound velocity, c , is related to the temperature
(iii) all the gas inside a sphere of accretion radius r S
centered at the cell is above a specified density thresholAdC,C T by the ideal gas law c2S =γRgasT/µ, where γ is the ef-
fective barotropic exponent, R is the ideal gas constant,
(iv) the gas inside the sphere is gravitationally bound, gas
and µ is the mean molecular weight. Inside the layer, we
Jeans-unstable and converging,
set γ = 1.0001 (effectively isothermal, γ = 1.0 is excluded
(v) radius r of the sink particle would not overlap
ACC with the adopted numerical scheme) µ=2 (pure molecular
with accretion radius of an already existing sink particle.
hydrogen), and T =10K.
Following(Federrathetal.2010),wesettheaccretionradius
to 2.5 grid cell size at the highest refinement level.
3.3 The external medium
All the simulations have mixed boundary conditions
(BCs)forself-gravity,i.e.periodicinthetwodirectionsx,y We use two different kinds of the layer confinement: in Sec-
parallel to the layer, and isolated in the third direction z tions4,5and6,weinvestigatelayersconfinedwiththermal
perpendicular to the layer. In order to calculate the gravi- pressurefrombothsurfaces;inSection7,weinvestigatelay-
tational field in this configuration, it is natural to seek for ers confined with the ram pressure from one surface, and
a modification ofthe standardEwaldmethod(Ewald1921; with thermal pressure from the other.
(cid:13)c 2017RAS,MNRAS000,1–21
6 F. Dinnbier et al.
The medium imposing the thermal pressure is imple- inducestwoartificialeffects:asmallrampressureactingon
mented as follows. In order to compare our simulations the contact discontinuity and an increase of layer surface
with the analytic theory, it is necessary that the ambient density Σ. We discuss the influence of the ram pressure at
medium has no dynamical influence on the layer, apart the end of Section 7.3 and show that it is not important.
from exerting the thermal pressure. By means of conver-
gence tests we establish that this can be achieved by set-
ting ρ (z ) (cid:46) 10−2ρ . In addition, to diminish the 4 LAYERS CONFINED BY THERMAL
AMB MAX O
influenceofsoundwavesreflectionfrombordersofthecom- PRESSURE WITH MONOCHROMATIC
putational domain, we extend the computational domain PERTURBATIONS
to ∼ ±3z . Therefore, we set the ambient medium to
MAX In this section, we test the dispersion relations by compar-
be isothermal with temperature T = 300K and mean
AMB ing them to simulations. We study two extreme cases of
molecular weight µ = 0.6. To prevent the ambient
AMB pressureconfinement:self–gravitydominatedwithA=0.99
medium from falling on the layer, we set its density profile
andexternalpressuredominatedwithA=0.18.Theapplied
close to the hydrostatic equilibrium, i.e.
perturbation is of a single wavelength (monochromatic).
ρ (z) = ρ (z )× Since the analytical estimates are based on linearised
AMB AMB MAX
(cid:18) 2πGΣ µ (|z|−z )(cid:19) equations,theyarevalidonlyaslongastheperturbingam-
exp − 0 AMB MAX . (13)
RgasTAMB plitude q1 of any quantity is smaller than its unperturbed
value q . Accordingly, we define the linear regime of frag-
In the case of accreting layers, ram pressure is realised O
mentationifthemaximumoftheperturbedsurfacedensity
by supersonic accretion of gas with uniform density ρ ,
temperatureT andsoundvelocityc .ThegasimpAaCcCts Σ1 issmallerthanΣO andthenon–linearregimeotherwise.
ACC ACC Thedispersionrelationcanbedeterminedonlyinthelinear
at velocity v (and hence Mach number M = v /c )
ACC ACC S part of the fragmenting process.
onto the upper surface of the layer at z (cid:39) z . The ac-
MAX The generic name of a monochromatic simulation is in
creted gas is cold and molecular and the shock is isother-
theformM<A> <kH >,wherethefirsttwonumbersaf-
mal, i.e. T = T and c = c . The bottom surface at HT
ACC ACC S ter ”M” represent the value of parameter A and the num-
z (cid:39) −z is confined by the ambient medium at temper-
MAX bersaftertheunderscoretheperturbingwavenumberinthe
ature T , T (cid:29) T exerting the thermal pressure on
AMB AMB dimensionless form kH . Thus, for example, simulation
the layer. To prevent the layer from bulk acceleration, we HT
M18 020treatsalayerwithA=0.18,andinitialmonochro-
set the pressures acting on both surfaces to be equal, i.e.
matic perturbation kH =0.20.
ρ (v2 +c2 )=ρ c2 . HT
ACC ACC ACC AMB AMB
The accretion leads to large scale flows inside the layer
(see Section 7.2 and Fig. 9), which mix the layer with the 4.1 Initial conditions for perturbations
warmerambientmedium.Consequently,ifcoolingwerenot
The initial monochromatic perturbation for a layer con-
implemented, the temperature of the layer would increase
fined by thermal pressure corresponds to the eigenfunction
and the layer would thicken. This is a spurious behaviour
for acoustic–surface–gravity modes of a thick layer with
which we suppress because a real layer would quickly cool
freely moving surfaces (see eqs. (20) – (23) in Kim et al.
to temperature T.
(2012)). The initial amplitude of perturbed surface density
To keep the layer at constant temperature, we need to
isΣ (0)=0.01Σ .Thelengthofthecomputationaldomain
distinguish it from the warm ambient medium. We detect 1 O
in direction x is equal to one perturbing wavelength.
thelayeraccordingtothesurfacedensityσ calculatedfrom
the bottom side z of the computational domain
BOT
4.2 The dispersion relation
z
(cid:90)
σ(z)= ρ(z(cid:48))dz(cid:48). (14) Simulations fortheexternalpressure dominated layer(A=
z 0.18,modelsM18)arelistedinTable1.Upperleftpanelof
BOT
Figure2showsevolutionofsurfacedensityperturbationsfor
We cool to temperature T any cell with the col-
selected modes. Since a small perturbing amplitude in the
umn density above threshold σCRIT. We set σCRIT = surface density behaves as Σ1(t)=Σ1(0)eiωt, the instanta-
−z
2 (cid:82)MAX ρ (z(cid:48))dz(cid:48)toenablethelayertofreelyripple.The neousgrowthrateω equalstotheslopeofthecurvesshown
AMB in the upper left panel of Fig. 2. ω is almost time indepen-
z
layeBrOdTetectionissensitivebecausethetotalcolumndensity dent throughout the simulations. Modes with kHHT<∼0.65
oftheambientmediumisoftheorderofthecolumndensity are unstable and grow, modes with kHHT>∼0.65 are stable.
ofonecellinsidethelayer,soσ(z)risessteeplyoncezenters Wemeasureω bylinearfittoln(Σ1/Σ1(0))overtimeinter-
the layer. val tfit, tfit. The starting time tfit is determined so as to
O 1 O
Althoughmostofthemassdeliveredtothelayercomes suppresstheinfluenceoftheinitialconditions(whichareex-
from the accreted medium, the intermixing consumes a sig- actly that of E78 eigenvectors). We suppose that the initial
nificant amount of the warm ambient medium. To prevent growth rate should be significantly altered at the timescale
the ambient medium from being exhausted in the course of of sound crossing time through one half of the wavelength,
a simulation, we continuously replenish gas at temperature i.e. (cid:39) πH /c (we use tfit = 0.1Myr). The upper bound
HT S O
T through the bottom boundary of the computational tfit is constrained by the condition for the perturbation to
AMB 1
domain,sothatthetotalmassofambientgasremainscon- be small, i.e. Σ < Σ . Since model M18 terminates be-
1 O
stant. The fresh ambient gas moves towards the layer and fore this condition is fulfilled, we set tfit near the end of
1
(cid:13)c 2017RAS,MNRAS000,1–21
Fragmentation of stratified gaseous layers 7
Run A nx×ny×nz HHT/dz kHHT tEEF78OLD tNUM
[Myr] [Myr]
M18 020 0.18 640×128×128 20.2 0.20 0.16 0.16
M18 025 0.18 512×128×128 20.2 0.25 0.15 0.16
M18 033 0.18 384×128×128 20.2 0.33 0.15 0.16
M18 050 0.18 256×128×128 20.2 0.50 0.19 0.19
M18 062 0.18 208×128×128 20.2 0.62 0.41 0.41
M18 073 0.18 176×128×128 20.2 0.73 - -
M18 089 0.18 144×128×128 20.2 0.89 - -
M99 025 0.99 512×128×128 20.2 0.25 0.74 0.76
M99 033 0.99 384×128×128 20.2 0.33 0.69 0.70
M99 050 0.99 256×128×128 20.2 0.50 0.67 0.68
M99 073 0.99 176×128×128 20.2 0.73 0.77 0.78
M99 089 0.99 144×128×128 20.2 0.89 1.12 1.14
M99 114 0.99 112×128×128 20.2 1.14 - -
M99 133 0.99 96×128×128 20.2 1.33 - -
Table 1. Parameters for simulations of layers confined by thermal pressure with monochromatic perturbations (models M). We list
parameterA,numberofcellsatthehighestrefinementlevelnx,ny andnz,resolutionintheverticaldirectionHHT/dz,parameterkHHT
where k is the selected wavenumber, tEFOLD is the analytical e–folding time according to E78 for given k, and t (t =1/ω) is
E78 NUM NUM
thee–foldingtimemeasuredinoursimulations.
Figure2.Thedispersionrelationforaself–gravitydominated(A=0.99;toprow)andexternal–pressuredominated(A=0.18;bottom
row) layer. Left panels: Time evolution of the surface density perturbations for monochromatic models. Only selected models are
plottedtoavoidconfusion.ThevalueofparameterkH foraparticularmodelisonrightfromthecurve.Thegrowthrateiscalculated
HT
intheintervalmarkedbytheverticaldottedlines.Rightpanels:Comparisonbetweenanalytical(E78,V83andW10)andnumerically
obtained dispersion relations. The growth rate obtained from monochromatic and polychromatic simulations is plotted by crosses and
dots,respectively.Dataforpolychromaticsimulationsisbinnedtoreducenoise.SeeSections4.2and5.2foradetaileddescription.
the simulation. We compare measured ω with the analyti- ations towards the end only when the perturbing ampli-
cal dispersion relations Eq. (7), Eq. (9) and Eq. (10) in the tude becomes large Σ >Σ . Measured ω (tfit = 1.4Myr,
1∼ O O
upper right panel of Figure 2 and list measured e–folding tfit = 3.0Myr) is in a good agreement with E78 disper-
1
time t = 1.0/ω alongside the E78 analytical estimate sionrelationandinconsistentwithV83(lowerrightpanelof
NUM
tEFOLD in Table 1. Our results are very close to E78 (solid Fig. 2). Since W10 and E78 are similar, our data could not
E78
line) and are inconsistent with both W10 (dotted line) and distinguishbetweenthemforself–gravitydominatedlayers.
V83 (dashed line).
Monochromatic simulations for the self–gravity domi-
nated layer (A = 0.99, models M99) also show nearly time
independent ω (lower left panel of Figure 2) with devi-
(cid:13)c 2017RAS,MNRAS000,1–21
8 F. Dinnbier et al.
5 LAYERS CONFINED BY THERMAL of Figure 2) as for the monochromatic models presented in
PRESSURE WITH POLYCHROMATIC Section4.Thedataarebinnedtoreducenoise.Thecritical
PERTURBATIONS wavenumber k as well as the growth rate of wavenum-
MAX
bersk<k areagaininaverygoodagreementwithE78
In this Section, we investigate fragmentation of self– MAX
both in self–gravity and external pressure dominated cases.
gravitatingpressureconfinedlayersbothinlinearandnon–
Wavenumbers with k > k are stable. We do not detect
linear regime. Initial perturbations contain many wave- MAX
any significant deviation between our results and E78.
lengthssimultaneously(polychromatic).Inthelinearregime
of fragmentation, we measure the dispersion relation. We
also study the fragmenting process qualitatively as a func- 5.3 Evolution in non–linear regime
tionofparameterA,andcomparethefragmentmassesand
Evolution of surface density for self–gravity dominated
fragmenting timescales with analytical estimates.
model P99 is shown in Figure 3. The fragmentation be-
The generic name of a polychromatic simulation con-
gins with emergence of round objects (at time from 0.7Myr
sists of letter ”P” followed by two numbers indicating the
to 2.1Myr). These objects then gradually grow and become
value of the parameter A. For example, simulation P18
slender as was predicted by the second order perturbation
treats a layer with A=0.18.
theorybyMiyamaetal.(1987b).Thetransformationofob-
jectsfromroundishtofilamentary–likeisapparentbetween
5.1 Initial conditions for perturbations plotsat2.8Myrto4.9Myr.Sinkparticlesforminthedensest
partsofthefilamentsandaccretematerialfromsurrounding
Toinitiatesimulationsoflayersconfinedbythermalpressure
filaments.Comparingthelastplotat4.9Myrwithaplotat
withpolychromaticperturbations,weperturbthelayerwith
the early stage of fragmentation (e.g. 2.1Myr), we see that
manywavevectorspointinginallthreespatialdirections.In
whenaclumpappears,itmonolithicallycollapsestoagrav-
order to be able to perform resolution tests, we generate
itationally unstable object.
their amplitudes in the Fourier space and then map them
Fragmentation in pressure–dominated case is repre-
bytheinverseFouriertransformonagrid,sotheirspectrum
sented by model P18 (Fig. 4). At the beginning of frag-
doesnotdependongridresolution.TheamplitudesA˜(k)are
mentation,thelayerswiftlybreaksintosmallobjects(plots
drawn as random variables from the uniform distribution
from0.3Myrto1.0Myr)withmassessmallerthantheJeans
for (cid:107)(cid:126)k(cid:107) < k and are zero otherwise, so the amplitudes
O mass M . The subjeans masses are a direct consequence
occupy a sphere in the Fourier space. To assure reasonable J
of fragmentation according to E78 for layers with low A
resolution, the highest wavenumber k corresponds to the
O (right panel of Fig. 1). For confinement of the objects, ex-
size of at least 4 grid cells. The whole unstable range of
ternalpressureismoreimportantthanself–gravity,soapart
wavenumbers predicted by any of the discussed estimate to
frombeingimmersedinanexternalgravitationalfieldofthe
the dispersion relation is always included inside the sphere,
layer, the objects are equivalents to gravitationally stable
i.e. max(k ,k ,k )<k .
E78.MAX V83.MAX W10.MAX O Bonnor–Ebert spheres. The stable objects then gradually
In order to compare fragmenting timescales between
mergeuntilenoughmassforagravitationallyboundclump
layerswithdifferentvaluesofparameterA,itisnecessaryto
is assembled (see plots from 1.0Myr to 2.3Myr). Merging
choose comparable amplitudes of their initial perturbation.
oftenleadstonon–radialaccretionresultinginspinning–up
To reach the aim, we normalise all the perturbing ampli-
(cid:113) thefragmentsanddiscformationaroundthem(plotattime
tudesA˜(k)tosatisfy (cid:104)A˜(k)2(cid:105)/A˜(0)=const.Tobeableto 2.3Myr). As a bound clump is formed, it collapses and its
calculate the dispersion relation in the simulations, the ini- cross section for possible following mergers is reduced and
tial perturbing amplitudes must be small. This imposes an merging rate decreases. Therefore the fragmenting process
upper limit on the normalisation constant. The lower limit in the pressure dominated case, which proceeds via coales-
is constrained by accuracy of the tree code. To fulfill these cence of many small clumps, is qualitatively different from
requirements, we set Σ1 (cid:39)0.1ΣO. the continuous collapse in the self–gravity dominated case.
We refer to the former and latter as coalescence driven col-
lapse and monolithic collapse, respectively.
5.2 The dispersion relation
Since the analytical dispersion relations are often used
To measure the dispersion relation, we compute Fourier for estimating fragmenting time–scales and mass of frag-
transformofsurfacedensityforanyframeofthesimulation ments, it is interesting to compare these quantities with
andthenfitgrowthrateofeachindividualmodeintimein- those found in our simulations. We identify gravitationally
terval(tfit,tfit).Astheinitialconditionsdonotcorrespond boundfragmentswiththealgorithmdescribedinAppendix
O 1
totheeigenfunctions,themodesrelaxatthebeginning,and A. We experiment with several definitions of fragmenting
need some time before start growing at a temporarily con- time t , and conclude that taking t to be an instant
FRG FRG
stant growth rate. The choice of tfit is constrained so as whenthetotalmassofgravitationallyboundobjectsexceeds
O
themodeswiththelongestwavelengthspresentinthecom- 1/2 of the total mass of the layer is a reasonable estimate
putational domain already grow at a constant rate. tfit is for following reasons. As formation of bound objects starts,
1
constrained so as the modes with the highest growth rate theirtotalmassincreasesrapidly,sotakinganotherfraction
have not reached the nonlinear regime yet. We find that of the total mass than 1/2 would not lead to a significantly
tfit =1/ω and tfit =3/ω fulfill at best the differenttimescale.Timet isalsoclosetothetimet
O E78.FAST 1 E78.FAST FRG SINK
requirements for our parameter choice. whensinkparticlesintotalcontain1/2ofthetotalmassof
Wecalculatethedispersionrelationforthesamevalues the layer (see Table 2 and left panel of Fig. 5).
ofparameterA(modelsP18andP99inTable2;rightpanels Comparisonbetweent andE78estimatetEFOLD =
FRG E78.FAST
(cid:13)c 2017RAS,MNRAS000,1–21
Fragmentation of stratified gaseous layers 9
Run A nx×ny×nz HHT/dz NJ MJ M¯CL tEEF78O.FLADST tFRG tSINK tNL
[M ] [M ] [Myr] [Myr] [Myr] [Myr]
(cid:12) (cid:12)
P18 0.18 1024×1024×64 4.0 8.8 2.3 1.1 0.15 2.01 2.15 0.42
P20 0.20 1024×1024×64 4.4 9.3 2.5 1.9 0.17 2.19 2.40 0.50
P22 0.22 1024×1024×64 4.4 12.8 2.8 1.7 0.18 2.11 2.23 0.58
P25 0.25 1024×1024×64 4.4 19.0 3.1 3.1 0.21 1.90 2.20 0.80
P30 0.30 512×512×64 3.9 9.8 3.6 2.1 0.25 2.20 2.40 0.99
P40 0.40 512×512×64 4.7 17.0 5.0 3.5 0.33 2.38 2.64 1.45
P60 0.60 512×512×64 6.6 29.7 7.4 10.0 0.47 3.00 3.35 2.22
P80 0.80 512×512×64 5.9 89.0 9.9 20.6 0.59 3.60 4.00 2.31
P99 0.99 512×512×64 9.1 71.0 12.5 32.0 0.66 4.15 4.90 2.90
Table 2. Parameters for simulations of layers confined by thermal pressure with polychromatic perturbations (models P). First four
columns have the same meaning as the columns in Table 1. Further, we provide number of Jeans masses in the computational domain
NJ,JeansmassMJ,meanmassofgravitationallyboundobjectsM¯CL atfragmentingtimetFRG,timetSINK,analyticale–foldingtime
tEFOLD forthemostunstablewavenumberk andtransitiontimebetweenlinearandnon–linearregimet .
E78.FAST FAST NL
Figure 3.Evolutionofsurfacedensityfortheself–gravitydominatedlayer(modelP99).Numberattheupperleftcorneristhetimein
Myr.Sinkparticlesareplottedaswhitecircles.Areaofaparticularcirclecorrespondstothesinkparticlemass.Notethatthecolourscale
forupperandbottomrowisdifferent.
1/ω (Eq.(11))asafunctionofparameterAisplot- e–folding times when fragmentation happens n =
E78.FAST EFOLD
tedintheleftpanelofFig.5.Wetestthecommonlyadopted t /tEFOLD isonlyafunctionofAforanylayer.Ifweob-
FRG E78.FAST
assumption that the fragmentation occurs at a constant tainn (A)fromoursimulations,wecanestimatet
EFOLD FRG
number of analytic e–folding times tEFOLD regardless of for any layer because tEFOLD is already known from Eq.
E78.FAST E78.FAST
A,sowemultiplytEFOLD byaconstanttomatchthedata (11).
E78.FAST
point for run P99. The assumption of constant number of
Thedependenceofn onAisshowninthemiddle
e–folding times holds in self–gravity dominated case, but EFOLD
panel of Fig. 5 (solid line). Whereas n is nearly con-
it underestimates the collapse time when external pressure EFOLD
stantintheself–gravitydominatedcase,itstronglyincreases
dominates.
as the external pressure increases. Before drawing conclu-
The parameters describing the layer are Σ , c and
O S sionsfromthisresult,weshouldverifythatitisnotsimply
P .Fromdimensionalanalysisitfollowsthatanyquantity
EXT caused by the fact that self–gravity dominated models al-
with the dimension of time is given by
readystartwitheffectivelyhigherinitialperturbations.For
c f(A)
t= S , (15) any model, the transition to the non–linear regime at time
GΣ
O tNL (when Σ1 = ΣO) occurs at almost constant number of
where f(A) is an unknown function of the parameter A. tEFOLD (Table2andmiddlepanelofFig.5),sotheinitial
E78.FAST
Since both tEFOLD and t are given by Eq. (15) (with perturbations are in this sense comparable among models
E78.FAST FRG
presumably different description of f(A)), the number of with different A. Taking t instead of t = 0 as the time
NL
(cid:13)c 2017RAS,MNRAS000,1–21
10 F. Dinnbier et al.
Figure 4.Evolutionofsurfacedensityfortheexternalpressuredominatedlayer(modelP18).CaptionisthesameasinFig.3.
whentheperturbingamplitudesarecomparablewouldeven thedecreaseoftheJeansmasswithincreasingexternalpres-
emphasizethedependenceofthenumberoftEFOLD onA; sure.
E78.FAST
when A is near unity, the fragmentation occurs soon after
t ,whenAissmall,ittakesmanye–foldingtimestoreach
NL
t . This behaviour reflects the two different fragmenting
FRG
scenarios;whenAisnearunity,theclumpscontinuecollaps-
ing, when A is small, the clumps are stable and gradually
merge, which causes the delay in fragmentation time. We
conclude that neither E78 describes fragmenting time cor- 6 LAYERS CONFINED BY THERMAL
rectly because fragmentation does not occur at a constant PRESSURE AND POSSIBLE PATTERN
number of e–folding times for various A. FORMATION IN SURFACE DENSITY
The mean mass of bound objects M¯ at the frag- Fuchs (1996) proposes that when fragmentation of an in-
CL
menting time as a function of A is listed in Table 2 and finitelythindiscbecomesnon–linear,atripleofmodeswith
shown in the right panel of Figure 5. Analytical estimates wavevectors(cid:107)(cid:126)k(cid:107)(cid:39)kFAST inclinedatanglesaround60◦ (i.e.
for V83 and E78 (M and M , respectively) and the intheFourierspaceforminganequilateraltriangle)hasthe
V83 E78
Jeans mass M for the midplane density (we use M = highest growth rate due to interaction. Wu¨nsch & Palouˇs
4.26 Ac4/(G2ΣJ ) from eq. (13-33) in Spitzer (1978))Jare (2001) semi–analytically find similar behaviour for the sur-
plotted bsy linesO. The fragment masses for E78 and V83 are face of a shell assuming V83. As a result, the modes create
estimatedfromthewavelengthwiththehighestgrowthrate a hexagonal pattern in surface density. On the other hand,
(i.e. M =πΣ (λ /2)2 where λ =2π/k ). In the applying second–order perturbation theory, Miyama et al.
self–gravity dOomiFnAaStTed case, the mFAaSsTs of fragmFAeSnTts is in a (1987a)showthatafragmentinglayerbreaksintogradually
goodagreementwiththeE78prediction.Thisresultisanat- slenderingfilamentswithnosignsofthehexagonalpattern.
ural consequence of the monolithic collapse since the mass Inthissection,wetesttheFuchs’propositionbysearch-
of fragments formed in the linear regime does not signifi- ingregularpatternsinsurfacedensityinourpolychromatic
cantly change later on, during the non–linear collapse. The simulations. Since we find no evidence for any pattern, we
thin–shellapproximationsystematicallyunderestimatesthe performthreededicatedmodelstotestwhetherpatternfor-
fragment masses for A>0.6. mation arises at all under very idealised circumstances.
∼
Name of the special models exploring possible pattern
In the external pressure dominated case, the mass of formation is in the form I<A> <N > [S|D], where the
w
fragments is of the order of the Jeans mass. It is a result first two numbers after ”I” represent the parameter A and
of their formation process, since the fragments are initially the number after underscore the number N of amplified
w
subjeansandstableuntiltheyassembledapproximatelythe wavepackets. For models with suffix S, the wavepackets are
Jeansmassandthencollapse.Therefore,theestimatebased amplified by a smooth function (wavepackets and the func-
on E78 leads to too small fragment mass. In contrast, V83 tion are described in Section 6.1), while for models with
overestimatesfragmentmassasitdoesnottakeintoaccount suffix D, the wavepackets are degenerated to single modes.
(cid:13)c 2017RAS,MNRAS000,1–21