TY - JOUR
AU - Palouš,, Jan
AB - Abstract 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 surface and ram pressure on the other. In the linear regime of fragmentation, the dispersion relation we obtain agrees well with that derived by Elmegreen & Elmegreen, and consequently deviates from the dispersion relations based on the thin shell approximation or pressure assisted gravitational instability. In the non-linear regime, the relative importance of the confining pressure to the self-gravity is a crucial parameter controlling the qualitative course of fragmentation. When confinement of the layer is dominated by external pressure, self-gravitating condensations are delivered by a two-stage process: first the layer fragments 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 small contribution to confinement of the layer, the layer fragments monolithically into gravitationally unstable clumps and there is no coalescence. This dichotomy persists whether the external pressure is thermal or ram. We apply these results to fragments forming in a shell swept up by an expanding H ii region, and find that, unless the swept-up gas is quite hot or the surrounding medium has low density, the fragments have low mass (⁠|${\lesssim} 3\, \mathrm{M}_{_{\odot}}$|⁠), and therefore they are unlikely to spawn stars that are sufficiently massive to promote sequential self-propagating star formation. hydrodynamics, instabilities, waves, stars: formation, H ii regions, ISM: kinematics and dynamics 1 INTRODUCTION Massive stars (⁠|$M_\star \gtrsim 8\, \mathrm{M}_{_{\odot }}$|⁠) strongly influence the interstellar medium (ISM) surrounding them, mainly via photoionization, stellar winds and supernova explosions. Elmegreen & Lada (1977) propose a mechanism (collect and collapse) whereby an overpressured H ii region, driven by young massive stars, expands into dense molecular gas. The expansion induces a spherical shock, and the surrounding gas accumulates between the shock and the ionization front. The resulting shell of cold gas increases in mass, and eventually fragments to form a new generation of stars. A crucial issue is the maximum mass of these newly formed stars. If some of them are sufficiently massive to excite a new H ii region, the process can repeat recursively and star formation propagates itself sequentially. Otherwise, star formation is quenched. Various triggering mechanisms are discussed in Elmegreen (1998). From an observational perspective, shells are common morphological structures in the ISM. Ehlerová & Palouš (2005) have detected more than 600 shells in H i. Churchwell et al. (2007) have identified 322 complete or partial rings in the infrared. Simpson et al. (2012) list more than 5000 infrared shells. Deharveng et al. (2010) find that at least 86 per cent of the infrared shells identified by Churchwell et al. (2007) encircle H ii regions ionized by O- and early B-type stars, suggesting that the shells are due to feedback from these stars. After a careful examination, Deharveng, Zavagno & Caplan (2005) find 17 infrared shells that are candidates for the collect-and-collapse mechanism. Evidence for propagating star formation has also been reported in significantly larger H i shells (Dawson et al. 2008, 2011; Egorov et al. 2014), indicating that feedback operates over a large range of scales. Further observational support for propagating star formation comes from age sequences of OB associations and star clusters in the vicinity of star-forming regions (e.g. Blaauw 1964, 1991; Brown, de Geus & de Zeeuw 1994; Bik et al. 2010) In order to estimate the properties of fragments condensing out of a swept-up shell, Elmegreen & Elmegreen (1978, hereafter E78), Doroshkevich (1980) and Lubow & Pringle (1993) investigate the stability of a vertically resolved self-gravitating layer,1 and derive a semi-analytical formula for the corresponding dispersion relation. Vishniac (1983, hereafter V83) derives the dispersion relation for a self-gravitating infinitesimally thin shell, and this has been used to estimate the properties of clumps condensing out of fragmenting shells. Whitworth et al. (1994a) and Elmegreen (1994) analyse the fragmentation of accreting shells, while Wünsch & Palouš (2001) consider non-accreting shells expanding into a vacuum. These dispersion relations are based on linear perturbation theory, so they are relevant only as long as the perturbing amplitudes are small. Miyama, Narita & Hayashi (1987a,b) extend the linearized solution of Elmegreen & Elmegreen (1978) into the non-linear regime, by including second-order terms. They find that a layer breaks into filaments that become increasingly slender with time. In contrast, Fuchs (1996) proposes that, in the non-linear regime, the fragments form a semi-regular hexagonal pattern on the surface of the layer. In order to test the validity of these analytic derivations, Dale et al. (2009) obtain a dispersion relation, based on numerical simulations of expanding, non-accreting shells, confined by constant external pressure. They find significant differences between their dispersion relation and that of V83 (based on the thin shell approximation). In subsequent work, Wünsch et al. (2010, hereafter W10) explain the difference by modifying V83 to include the effect of pressure confinement. They call the mechanism underlying their dispersion relation pressure assisted gravitational instability (PAGI; Wünsch et al. 2010). Since Dale et al. (2009) simulate fragmentation of a whole shell, and the shell gets thinner with increasing external pressure, they are only able to resolve shells that are confined by low to moderate values of the external pressure, and therefore the W10 dispersion relation is unverified in the case of high external pressure. However, it is in the limit of high external pressure that W10 differs substantially from E78. Van Loo, Keto & Zhang (2014) investigate self-gravitating layers permeated by magnetic fields and present a non-magnetic control run. Although their results differ from V83, they are in agreement with both E78 and W10, because their layer is confined by a very low external pressure. Iwasaki, Inutsuka & Tsuribe (2011b) also derive a numerical dispersion relation that differs significantly from that of V83. Our main aims are to simulate the fragmentation of gaseous layers, in both the low and high ambient pressure cases, and to compare the results with the analytic or semi-analytic estimates derived in previous studies. This comparison addresses three issues: (i) dispersion relations in the linear regime of fragmentation; (ii) the elapsed time before the layer forms gravitationally bound fragments, and the resulting fragment masses; and (iii) the possibility that mode interaction leads to fragments distributed on a regular periodically repeating pattern. We study only a small square patch on the layer, so that we can attain good resolution in directions perpendicular to the layer, even when the layer is significantly compressed. In addition to layers confined from both sides by thermal pressure, we also simulate layers accreting on to one surface. We use these results to investigate the fragmentation of a shell driven by an expanding H ii region – using an analytic solution to account for the expansion of the H ii region – and contrast our results with those obtained analytically by Whitworth et al. (1994a), and by Iwasaki et al. (2011b) who performed simulations with small perturbing amplitudes. The paper is organized as follows. Section 2 reviews the most important properties of self-gravitating layers and the analytical estimates (E78, V83, W10) of the dispersion relations. These estimates are used for comparison with the numerical results. Section 3 describes applied numerical methods and initial conditions. Section 4 describes simulations of layers confined on both sides by thermal pressure (due to a very hot rarefied gas), and seeded with monochromatic perturbations. Section 5 describes simulations of layers confined on both sides by thermal pressure, and seeded with polychromatic perturbations. Section 6 describes simulations in which we explore whether fragments forming in layers tend to be arranged into regular patterns. Section 7 describes simulations of layers in which one side is confined by thermal pressure and the other by the ram pressure of a homogeneous plane-parallel inflow. Section 8 considers an H ii region expanding into a homogeneous medium, and estimates the time at which the swept-up shell fragments, and the properties of the fragments. Section 9 discusses the results, and Section 10 summarizes our main conclusions. Appendix A describes our algorithm developed for finding gravitationally bound fragments. 2 LINEARIZED THEORY OF LAYER FRAGMENTATION 2.1 The unperturbed state Before introducing the dispersion relations describing the growth rate of perturbations in the self-gravitating pressure-confined layer, we review the unperturbed configuration. The model is also used to generate the initial conditions for our simulations. We assume that the layer is initially in plane-parallel stratified hydrostatic equilibrium, i.e. it is infinite in the x and y directions, its normal points in the z direction and all quantities are functions of z only. The layer is isothermal with sound speed cS, so the density distribution is as derived by Spitzer (1942) and Goldreich & Lynden-Bell (1965), i.e. \begin{equation} \rho (z)=\frac{\rho {_{\rm O}}}{\cosh ^2(z/H{_{\rm O}})}, \end{equation} (1) where ρO is the density on the mid-plane (z = 0), \begin{equation} H{_{\rm O}}=\frac{c{_{\rm S}}}{\sqrt{2 \pi G \rho {_{\rm O}}}} \end{equation} (2) is the vertical scaleheight and G is the gravitational constant. The surfaces of the layer are at z = ±zMAX, and outside this (|z| > zMAX) there is a hot gas with negligible density that exerts an external pressure PEXT. The layer is fully characterized by PEXT, cS and its surface density ΣO. A dimensionless parameter constructed from these quantities (Elmegreen & Elmegreen 1978), \begin{equation} A=\frac{1}{\sqrt{1+\left(2 P{_{\rm EXT}}/\pi G \Sigma {_{\rm O}}^2\right)}}\,, \end{equation} (3) reflects the relative importance of self-gravity and external pressure in holding the layer together. Parameter A increases monotonically from nearly 0 (external pressure dominated; PEXT ≫ GΣO2) to 1 (self-gravity dominated). From pressure equilibrium at the surfaces, PEXT = ρ(zMAX)cS2, it follows \begin{equation} z{_{\rm MAX}}=H{_{\rm O}}\mathrm{arcosh}^{-1}\!\left( \frac{\rho {_{\rm O}}c{_{\rm S}}^2 }{P{_{\rm EXT}}}\right)=H{_{\rm O}}\mathrm{arcosh}^{-1}\!\left( \frac{1}{\sqrt{1-A^2}}\right)\!.\!\!\!\!\!\! \end{equation} (4) Half-thickness of a stratified layer is defined by HHT ≡ ΣO/(2ρO). Equations (1)–(3) then yield \begin{equation} H{_{\rm HT}}= \frac{\Sigma {_{\rm O}}c{_{\rm S}}^2}{2P{_{\rm EXT}}+ \pi G \Sigma {_{\rm O}}^2} = \frac{c{_{\rm S}}^2 A^2}{\pi G \Sigma {_{\rm O}}}\,, \end{equation} (5) and the mid-plane density, ρO, depends on ΣO as \begin{equation} \rho {_{\rm O}}=\frac{2P{_{\rm EXT}}+ \pi G \Sigma {_{\rm O}}^2}{2 c{_{\rm S}}^2} = \frac{\pi G \Sigma {_{\rm O}}^2}{2 c{_{\rm S}}^2 A^2}. \end{equation} (6) Pressure-dominated layers are of almost uniform density (zMAX ≃ HHT), while self-gravity-dominated layers have a pronounced density maximum on the mid-plane (zMAX ≫ HHT). 2.2 Analytical estimates of the dispersion relation In this section, we review and compare the assumptions underlying the dispersion relations derived by E78, V83 and W10. To simplify the discussion, we neglect the effects of shell curvature and velocity divergence, by setting the shell radius to infinity. The results are therefore applicable to a plane-parallel layer, and can also be applied to a shell, as long as the unstable wavelengths are much shorter than the radius of the shell. The dispersion relation gives the perturbation growth rate |$\omega$| as a function of wavenumber k. In planar geometry, and for the dispersion relations in question, |$\omega$| is either purely real or purely imaginary. In this paper, we adopt the convention that growing instability is described by the positive real part of |$\omega$|⁠, which we denote for simplicity ω, and we omit the oscillating imaginary part. The characteristic time-scale of perturbation growth, its e-folding time, is defined as tEFOLD = 1/ω. For all the dispersion relations in question, there is a finite range of unstable wavenumbers (0, kMAX), where kMAX is the highest unstable wavenumber. The maximum growth rate ωFAST is attained for wavenumber kFAST that is always approximately kMAX/2. The wavelengths corresponding to kMAX and kFAST are denoted by λMAX and λFAST, respectively. When describing a particular analytical estimate of the dispersion relation, the subscript begins with its name, e.g. kE78.FAST is the wavelength kFAST for E78. The abbreviations E78, V83 and W10 refer to the particular dispersion relation, not to the paper where they are first described. The unstable branches of the dispersion relations for both self-gravity- (A = 0.99) and external pressure-dominated (A = 0.18) layers are shown in the right column of Fig. 2. 2.2.1 Dispersion relation for the thin shell approximation To estimate the dispersion relation of an expanding shell, Vishniac (1983) reduce the problem to two dimensions by integrating the continuity, Euler and Poisson equation through the thickness of the shell. In planar geometry, the thin shell dispersion relation becomes \begin{equation} \omega_{{\rm V83}}^2(k) = 2\pi G \Sigma {_{\rm O}}k - c_{{\rm S}}^2 k^2. \end{equation} (7) The stability of modes is determined by the imbalance between self-gravity and internal pressure gradient, and there is no contribution from external pressure. 2.2.2 Dispersion relation for pressure assisted gravitational instability In order to evaluate the influence of external pressure, Wünsch et al. (2010) investigate perturbations with the form of an oblate spheroid, embedded in a layer. The spheroid is homogeneous and confined by external pressure PEXT. Its semi-major axis is r, its semi-minor axis is the layer's half-thickness, HHT, and its total mass is M. Radial excursion of an element on the equator of the spheroid is regulated by the equation of motion \begin{eqnarray} \nonumber \ddot{r} & = & -\frac{3GM}{2r^2} \left\lbrace \frac{\cos ^{-1}(H{_{\rm HT}}/r)}{(1-(H{_{\rm HT}}/r)^2)^{3/2}} -\frac{H{_{\rm HT}}/r}{1-(H{_{\rm HT}}/r)^2} \right\rbrace \\ && -\frac{20 \pi r H{_{\rm HT}}P{_{\rm EXT}}}{3 M} + \frac{5 c_{{\rm S}}^2}{r} \end{eqnarray} (8) (Boyd & Whitworth 2005). Wünsch et al. (2010) equate the instability growth rate to the radial contraction rate of the spheroid. Collapse from radius rO, by a small factor ε to radius (1 − ε)rO, in time tε, fulfils |$(\epsilon - 1)r{_{\rm O}}= \dot{r}{_{\rm O}}t_{\epsilon } + \frac{1}{2} \ddot{r}{_{\rm O}}t_{\epsilon }^2$|⁠. Using |$\omega =1/t_{\epsilon }$| and r = λ/2 = π/k, this yields \begin{eqnarray} \omega _{{\rm W10}}^2(k)&=&-\frac{\ddot{r}{_{\rm O}}}{2 \epsilon r{_{\rm O}}} \nonumber \\ &=&\frac{1}{\epsilon } \bigg \lbrace - \frac{5 c_{{\rm S}}^2 k^2}{2 \pi ^2} \nonumber \\ && +\, \frac{3\,G \Sigma {_{\rm O}}k}{4} \!\left(\!\frac{\cos ^{-1}(kH_{{\rm HT}}/\pi )\!-\!\sqrt{k^2 H_{{\rm HT}}^2 /\pi ^2\!-\!1}}{(1-k^2 H_{{\rm HT}}^2 /\pi ^2)^{3/2}}\right) \nonumber \\ && +\, \frac{10 P{_{\rm EXT}}c_{{\rm S}}^2 k^2}{3 \pi ^2 (2 P{_{\rm EXT}}+\pi G \Sigma _{{\rm O}}^2)} \bigg \rbrace . \end{eqnarray} (9) The terms on the right-hand side of equation (9) represent the gradient of internal pressure (second line), self-gravity (third line) and confinement by external pressure (fourth line). Wünsch et al. (2010) suggest setting ε ∼ 0.1, but, although ε affects the magnitude of |$\omega {_{\rm W10}}$|⁠, it does not influence the range of unstable wavenumbers. 2.2.3 Dispersion relation for a vertically stratified layer Elmegreen & Elmegreen (1978) obtain the dispersion relation for a self-gravitating, pressure-confined, semi-infinite layer in hydrostatic equilibrium by solving the system of perturbed continuity, Euler and Poisson equations. Kim et al. (2012) revisit their study and offer an additional insight. According to Kim et al. (2012), the dispersion relation can be written as \begin{equation} \omega _{{\rm E78}}^2 = 2 \pi G \Sigma {_{\rm O}}k (\eta F_{_{\rm J}}+(1-\eta )F_{_{\rm D}}) - c_{_{\rm EFF}}^2 k^2 \,. \end{equation} (10) Here |$c_{_{\rm EFF}}$| is the effective sound speed. |$F_{_{\rm J}}$| and |$F_{_{\rm D}}$| are reduction factors for self-gravity. Parameter η is the fraction of the perturbed surface density that is attributable to compression of material near the centre of a proto-fragment (and hence near the mid-plane of the layer). The rest of the perturbed surface density is due to corrugations on the surface of the layer. Thus, η controls the relative importance of compressional and surface-gravity waves. The former are important in self-gravity-dominated layers, while the latter are important in pressure-dominated layers. |$c_{_{\rm EFF}}$|⁠, η, |$F_{_{\rm J}}$| and |$F_{_{\rm D}}$| are complicated functions of k, which cannot generally be expressed in closed form. The growth rate ωE78.FAST of the most unstable wavenumber kE78.FAST in the limit A → 0 (layers dominated by external pressure) is (equation 47 in Kim et al. 2012) \begin{equation} \omega _{{\rm E78.FAST}}^2 = 0.276 \times 2 \pi G \rho {_{\rm O}}= 0.276 \left( \frac{\pi G \Sigma {_{\rm O}}}{c{_{\rm S}}A} \right)^2. \end{equation} (11) 2.2.4 Comparison between the analytical estimates of the dispersion relation The left-hand panel of Fig. 1 shows the dependence of kMAX normalized to kV83.MAX on parameter A for different dispersion relations. For self-gravity-dominated layers, both W10 and E78 predict almost identical ranges of unstable wavenumbers (kE78.MAX/kW10.MAX → 0.966 as A → 1), but V83 predicts a broader range with kW10.MAX/kV83.MAX → 0.518. With increasing external pressure (decreasing A), the difference between E78 and V83 diminishes. However, when A ≲ 0.5, E78 extends to much higher wavenumbers, indicating that shorter wavelengths are unstable. And, as A → 0, W10 and E78 have different limiting behaviour: kW10.MAX tends to a constant (kW10.MAX/kV83.MAX → 2.216), whereas kE78.MAX ∝ HHT−1 ∝ A−2. Figure 1. Open in new tabDownload slide Comparison between the dispersion relations derived by E78, V83 and W10, as a function of A.  Left-hand panel: the marginally stable wavenumber kMAX, for E78 and W10, normalized to kV83.MAX. Right-hand panel: mass of the fragment with the highest growth rate, MFAST, normalized to the mid-plane Jeans mass, MJ. Fragments formed with mass below ∼MJ (thin line) are gravitationally stable. Figure 1. Open in new tabDownload slide Comparison between the dispersion relations derived by E78, V83 and W10, as a function of A.  Left-hand panel: the marginally stable wavenumber kMAX, for E78 and W10, normalized to kV83.MAX. Right-hand panel: mass of the fragment with the highest growth rate, MFAST, normalized to the mid-plane Jeans mass, MJ. Fragments formed with mass below ∼MJ (thin line) are gravitationally stable. Since the wavelength λFAST has the highest growth rate, we assume that the fragment mass MFAST is approximately the mass confined inside a circle of radius λFAST/2, i.e. MFAST = πΣO(λFAST/2)2 = πΣO(π/kFAST)2. Masses MFAST normalized to the mid-plane Jeans mass MJ are plotted in the right-hand panel of Fig. 1. V83 and W10 predict that when formed, the fragments are already gravitationally unstable (MFAST/MJ ≳ 1) for any A. On the other hand, E78 predicts gravitationally unstable fragments only when A ≳ 0.5. According to E78, the layer breaks into gravitationally stable fragments for lower values of A, thus predicting qualitatively different scenario than V83 and W10. We simulate and discuss evolution of pressure-dominated layers (low A) in Sections 5.3 and 9.1. 3 METHODS AND INITIAL CONDITIONS 3.1 Numerics All the hydrodynamic simulations presented in the paper are performed with the MPI-parallelized flash4.0 code (Fryxell et al. 2000). flash is an AMR code based on the PARAMESH library (MacNeice et al. 2000). The hydrodynamic equations are solved by the piecewise parabolic method (Colella & Woodward 1984). Self-gravity is calculated using an octal tree code (Wünsch et al., in preparation) that offers three acceptance criteria for interaction between the target point (where the acceleration is evaluated) and a node (the source of the gravitational force). The algorithm invented by Barnes & Hut (1986), which accepts nodes seen from the target point at an angle smaller than a specified value. This algorithm is purely geometric since it does not take into account distribution of mass inside nodes, or their relative contribution to the net gravitational force. A set of criteria using the node size, mass and optionally higher multipole moments of the mass distribution within the node. They either estimate the upper limit on the acceleration error of the cell–node interaction Δamax from equation (9) of Salmon & Warren (1994) or they use the approximation by Springel (2005), \begin{equation} \Delta a_\mathrm{(p)}^\mathrm{max} = \frac{GM}{d^2}\left(\frac{h}{d}\right)^{p+1}, \end{equation} (12) where M and h are the node mass and size, respectively, d is the distance between the cell and the node mass centre, and p is the order of the multipole expansion. An experimental implementation of the ‘sumsquare’ criterion of Salmon & Warren (1994), which controls the sum of all errors’ origination from all individual node contributions. We invoke the criterion (ii) with equation (12) and p = 1 in all our simulations because our tests show that it provides the best compromise between the code performance and accuracy. We convert dense gaseous condensations into sink particles according to conditions described in Federrath et al. (2010). To create a sink particle inside a particular cell, all following conditions must be fulfilled: the cell is at the highest refinement level, the cell contains minimum of gravitational potential, all the gas inside a sphere of accretion radius rACC centred at the cell is above a specified density threshold, the gas inside the sphere is gravitationally bound, Jeans-unstable and converging, radius rACC of the sink particle would not overlap with accretion radius of an already existing sink particle. Following Federrath et al. (2010), we set the accretion radius to 2.5 grid cell size at the highest refinement level. All the simulations have mixed boundary conditions (BCs) for self-gravity, i.e. periodic in the two directions x, y parallel to the layer, and isolated in the third direction z perpendicular to the layer. In order to calculate the gravitational field in this configuration, it is natural to seek for a modification of the standard Ewald method (Ewald 1921; Klessen 1997). By computing the appropriate limit of the standard Ewald method, we find formulae for the gravitational acceleration and potential in closed form for a configuration with mixed BCs; the derivation is described in Wünsch et al. (in preparation). The modification is used to calculate the gravitational field in the simulations presented here. The hydrodynamic BCs are periodic in the x and y directions, and reflecting in the z direction for all but accreting simulations. For accreting simulations, the BCs are inflow from the top of the computational domain, and diode from the bottom to prevent reflections of waves. The column density needed for cooling the warm ambient gas intermixed with the cold layer during accreting runs (see Section 3.3 for details) is calculated using module TreeRay/OpticalDepth of the tree code. The grid cells are cubic in all the simulations. The half-thickness of the layer (HHT; see equation 5 below) is always at least four grid cells (see Tables 1–5). Consequently, the Jeans length is always resolved by more than four grid cells, and the simulations satisfy the criterion for avoiding artificial gravitational fragmentation, as given by Truelove et al. (1997). We have performed successful convergence tests throughout the range of conditions simulated. We have also checked that the portion of the layer inside the computational domain is sufficiently large, i.e. that the periodic copies in the x and y directions do not significantly influence the properties of fragments. 3.2 Initial conditions for the layer Each model starts with a layer with properties described in Section 2.1 seeded with a perturbation. The perturbation is either a single mode given by the eigenfunction for acoustic-surface-gravity modes or white noise. The exact form of the perturbation is described at the beginning of a corresponding section (Sections 4 to 7). Particular values for the layer parameters are adopted from Iwasaki et al. (2011b) who investigate an H ii region excited by a 41 M⊙ star, expanding into a medium with number density n = 103 cm−3. At time t = 0.81 Myr, the shell has radius R = 3.86 pc and surface density ΣO = 0.0068 g cm−2. We use this value of ΣO for the initial conditions of all the simulations presented here, and vary A by changing PEXT. However, we note that the model is sufficiently simple that its physical parameters, (ΣO, cS, PEXT), can be rescaled arbitrarily, as long as A is unchanged, i.e. PEXT ∝ ΣO2. The sound velocity, cS, is related to the temperature T by the ideal gas law cS2 = γRgasT/μ, where γ is the effective barotropic exponent, Rgas is the ideal gas constant and μ is the mean molecular weight. Inside the layer, we set γ = 1.0001 (effectively isothermal, γ = 1.0 is excluded with the adopted numerical scheme), μ = 2 (pure molecular hydrogen) and T = 10 K. 3.3 The external medium We use two different kinds of the layer confinement: in Sections 4–6, we investigate layers confined with thermal pressure from both surfaces; in Section 7, we investigate layers confined with the ram pressure from one surface and with thermal pressure from the other. The medium imposing the thermal pressure is implemented as follows. In order to compare our simulations with the analytic theory, it is necessary that the ambient medium has no dynamical influence on the layer, apart from exerting the thermal pressure. By means of convergence tests, we establish that this can be achieved by setting ρAMB(zMAX) ≲ 10−2ρO. In addition, to diminish the influence of sound wave reflection from borders of the computational domain, we extend the computational domain to ∼± 3zMAX. Therefore, we set the ambient medium to be isothermal with temperature TAMB = 300 K and mean molecular weight μAMB = 0.6. To prevent the ambient medium from falling on the layer, we set its density profile close to the hydrostatic equilibrium, i.e. \begin{eqnarray} \rho {_{\rm AMB}}(z) & = & \rho {_{\rm AMB}}(z{_{\rm MAX}}) \nonumber \\ && \times \, \exp \left(-\frac{2 \pi G \Sigma _0 \mu {_{\rm AMB}}(|z| - z{_{\rm MAX}})}{R_{{\rm gas}} T{_{\rm AMB}}} \right). \end{eqnarray} (13) In the case of accreting layers, ram pressure is realized by supersonic accretion of gas with uniform density ρACC, temperature TACC and sound velocity cACC. The gas impacts at velocity vACC (and hence Mach number |$\scr {M} = v{_{\rm ACC}}/c{_{\rm S}}$|⁠) on to the upper surface of the layer at z ≃ zMAX. The accreted gas is cold and molecular and the shock is isothermal, i.e. TACC = T and cACC = cS. The bottom surface at z ≃ −zMAX is confined by the ambient medium at temperature TAMB, TAMB ≫ T exerting the thermal pressure on the layer. To prevent the layer from bulk acceleration, we set the pressures acting on both surfaces to be equal, i.e. ρACC(⁠|$v_{\rm ACC}^2$| + |$c_{\rm ACC}^2$|⁠) = ρAMB|$c_{\rm AMB}^2$|⁠. The accretion leads to large-scale flows inside the layer (see Section 7.2 and Fig. 9), which mix the layer with the warmer ambient medium. Consequently, if cooling were not implemented, the temperature of the layer would increase and the layer would thicken. This is a spurious behaviour that we suppress because a real layer would quickly cool to temperature T. To keep the layer at constant temperature, we need to distinguish it from the warm ambient medium. We detect the layer according to the surface density σ calculated from the bottom side zBOT of the computational domain \begin{equation} \sigma (z)=\int _{z_{\rm BOT}}^{z}\,\rho (z^{\prime })\,{\rm d}z^{\prime }\,. \end{equation} (14) We cool to temperature T any cell with the column density above threshold σCRIT. We set |$\sigma {_{\rm CRIT}}= 2 \int _{z_{{\rm BOT}}}^{-z_{{\rm MAX}}}\,\rho {_{\rm AMB}}(z^{\prime })\,{\rm d}z^{\prime }$| to enable the layer to freely ripple. The layer detection is sensitive because the total column density of the ambient medium is of the order of the column density of one cell inside the layer, so σ(z) rises steeply once z enters the layer. Although most of the mass delivered to the layer comes from the accreted medium, the intermixing consumes a significant amount of the warm ambient medium. To prevent the ambient medium from being exhausted in the course of a simulation, we continuously replenish gas at temperature TAMB through the bottom boundary of the computational domain, so that the total mass of ambient gas remains constant. The fresh ambient gas moves towards the layer and induces two artificial effects: a small ram pressure acting on the contact discontinuity and an increase of layer surface density Σ. We discuss the influence of the ram pressure at the end of Section 7.3 and show that it is not important. 4 LAYERS CONFINED BY THERMAL PRESSURE WITH MONOCHROMATIC PERTURBATIONS In this section, we test the dispersion relations by comparing them to simulations. We study two extreme cases of pressure confinement: self-gravity dominated with A = 0.99 and external pressure dominated with A = 0.18. The applied perturbation is of a single wavelength (monochromatic). Since the analytical estimates are based on linearized equations, they are valid only as long as the perturbing amplitude q1 of any quantity is smaller than its unperturbed value qO. Accordingly, we define the linear regime of fragmentation if the maximum of the perturbed surface density Σ1 is smaller than ΣO and the non-linear regime otherwise. The dispersion relation can be determined only in the linear part of the fragmenting process. The generic name of a monochromatic simulation is in the form |$M{\langle }A\rangle {}\_{}{\langle }kH{_{\rm HT}}{\rangle }$|⁠, where the first two numbers after ‘M’ represent the value of parameter A and the numbers after the underscore the perturbing wavenumber in the dimensionless form kHHT. Thus, for example, simulation M18_020 treats a layer with A = 0.18, and initial monochromatic perturbation kHHT = 0.20. 4.1 Initial conditions for perturbations The initial monochromatic perturbation for a layer confined by thermal pressure corresponds to the eigenfunction for acoustic-surface-gravity modes of a thick layer with freely moving surfaces (see equations 20–23 in Kim et al. 2012). The initial amplitude of perturbed surface density is Σ1(0) = 0.01ΣO. The length of the computational domain in direction x is equal to one perturbing wavelength. 4.2 The dispersion relation Simulations for the external pressure-dominated layer (A = 0.18, models M18) are listed in Table 1. The upper-left panel of Fig. 2 shows evolution of surface density perturbations for selected modes. Since a small perturbing amplitude in the surface density behaves as Σ1(t) = Σ1(0)eiωt, the instantaneous growth rate ω equals to the slope of the curves shown in the upper-left panel of Fig. 2. ω is almost time independent throughout the simulations. Modes with kHHT ≲ 0.65 are unstable and grow, and modes with kHHT ≳ 0.65 are stable. We measure ω by linear fit to ln (Σ1/Σ1(0)) over time interval |$t^{\rm fit}_{\rm O}$|⁠, |$t^{{\rm fit}}_1$|⁠. The starting time |$t^{\rm fit}_{\rm O}$| is determined so as to suppress the influence of the initial conditions (which are exactly that of E78 eigenvectors). We suppose that the initial growth rate should be significantly altered at the time-scale of sound crossing time through one half of the wavelength, i.e. ≃ πHHT/cS (we use |$t^{\rm fit}_{\rm O}$| = 0.1 Myr). The upper bound |$t^{{\rm fit}}_1$| is constrained by the condition for the perturbation to be small, i.e. Σ1 < ΣO. Since model M18 terminates before this condition is fulfilled, we set |$t^{{\rm fit}}_1$| near the end of the simulation. We compare measured ω with the analytical dispersion relations (equations 7, 9 and 10) in the upper-right panel of Fig. 2 and list measured e-folding time tNUM = 1.0/ω alongside the E78 analytical estimate |$t_{{\rm E78}}^{{\rm EFOLD}}$| in Table 1. Our results are very close to E78 (solid line) and are inconsistent with both W10 (dotted line) and V83 (dashed line). Figure 2. Open in new tabDownload slide The dispersion relation for an external pressure-dominated (A = 0.18; top row) and self-gravity-dominated (A = 0.99; bottom row) layer. Left-hand panels: time evolution of the surface density perturbations for monochromatic models. Only selected models are plotted to avoid confusion. The value of parameter kHHT for a particular model is on right from the curve. The growth rate is calculated in the interval marked by the vertical dotted lines. Right-hand panels: comparison between analytical (E78, V83 and W10) and numerically obtained dispersion relations. The growth rate obtained from monochromatic and polychromatic simulations is plotted by crosses and dots, respectively. Data for polychromatic simulations are binned to reduce noise. See Sections 4.2 and 5.2 for a detailed description. Figure 2. Open in new tabDownload slide The dispersion relation for an external pressure-dominated (A = 0.18; top row) and self-gravity-dominated (A = 0.99; bottom row) layer. Left-hand panels: time evolution of the surface density perturbations for monochromatic models. Only selected models are plotted to avoid confusion. The value of parameter kHHT for a particular model is on right from the curve. The growth rate is calculated in the interval marked by the vertical dotted lines. Right-hand panels: comparison between analytical (E78, V83 and W10) and numerically obtained dispersion relations. The growth rate obtained from monochromatic and polychromatic simulations is plotted by crosses and dots, respectively. Data for polychromatic simulations are binned to reduce noise. See Sections 4.2 and 5.2 for a detailed description. Table 1. Parameters for simulations of layers confined by thermal pressure with monochromatic perturbations (models M). We list parameter A, number of cells at the highest refinement level nx, ny and nz, resolution in the vertical direction HHT/dz, parameter kHHT, where k is the selected wavenumber, |$t_{{\rm E78}}^{_{\rm EFOLD}}$| is the analytical e-folding time according to E78 for given k and tNUM (tNUM = 1/ω) is the e-folding time measured in our simulations. Run . A . nx × ny × nz . HHT/dz . kHHT . |$t_{{\rm E78}}^{{\rm EFOLD}}$| . 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 – – Run . A . nx × ny × nz . HHT/dz . kHHT . |$t_{{\rm E78}}^{{\rm EFOLD}}$| . 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 – – Open in new tab Table 1. Parameters for simulations of layers confined by thermal pressure with monochromatic perturbations (models M). We list parameter A, number of cells at the highest refinement level nx, ny and nz, resolution in the vertical direction HHT/dz, parameter kHHT, where k is the selected wavenumber, |$t_{{\rm E78}}^{_{\rm EFOLD}}$| is the analytical e-folding time according to E78 for given k and tNUM (tNUM = 1/ω) is the e-folding time measured in our simulations. Run . A . nx × ny × nz . HHT/dz . kHHT . |$t_{{\rm E78}}^{{\rm EFOLD}}$| . 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 – – Run . A . nx × ny × nz . HHT/dz . kHHT . |$t_{{\rm E78}}^{{\rm EFOLD}}$| . 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 – – Open in new tab Monochromatic simulations for the self-gravity-dominated layer (A = 0.99, models M99) also show nearly time independent ω (lower-left panel of Fig. 2) with deviations towards the end only when the perturbing amplitude becomes large Σ1 ≳ ΣO. Measured ω (⁠|$t^{\rm fit}_{\rm O}$| = 1.4 Myr, |$t^{{\rm fit}}_1 = 3.0 \, \mathrm{Myr}$|⁠) is in a good agreement with E78 dispersion relation and inconsistent with V83 (lower-right panel of Fig. 2). Since W10 and E78 are similar, our data could not distinguish between them for self-gravity-dominated layers. 5 LAYERS CONFINED BY THERMAL PRESSURE WITH POLYCHROMATIC PERTURBATIONS In this section, we investigate fragmentation of self-gravitating pressure-confined layers both in linear and non-linear regime. Initial perturbations contain many wavelengths simultaneously (polychromatic). In the linear regime of fragmentation, we measure the dispersion relation. We also study the fragmenting process qualitatively as a function of parameter A, and compare the fragment masses and fragmenting time-scales with analytical estimates. The generic name of a polychromatic simulation consists of letter ‘P’ followed by two numbers indicating the value of the parameter A. For example, simulation P18 treats a layer with A = 0.18. 5.1 Initial conditions for perturbations To initiate simulations of layers confined by thermal pressure with polychromatic perturbations, we perturb the layer with many wavevectors pointing in all three spatial directions. In order to be able to perform resolution tests, we generate their amplitudes in the Fourier space and then map them by the inverse Fourier transform on a grid, so their spectrum does not depend on grid resolution. The amplitudes |$\tilde{A}(k)$| are drawn as random variables from the uniform distribution for |$\Vert \boldsymbol {k}\Vert < k{_{\rm O}}$| and are zero otherwise, so the amplitudes occupy a sphere in the Fourier space. To assure reasonable resolution, the highest wavenumber kO corresponds to the size of at least four grid cells. The whole unstable range of wavenumbers predicted by any of the discussed estimate to the dispersion relation is always included inside the sphere, i.e. max(kE78.MAX, kV83.MAX, kW10.MAX) < kO. In order to compare fragmenting time-scales between layers with different values of parameter A, it is necessary to choose comparable amplitudes of their initial perturbation. To reach the aim, we normalize all the perturbing amplitudes |$\tilde{A}(k)$| to satisfy |$\sqrt{\langle \tilde{A}(k)^2 \rangle }/\tilde{A}(0) = {\rm const}$|⁠. To be able to calculate the dispersion relation in the simulations, the initial perturbing amplitudes must be small. This imposes an upper limit on the normalization constant. The lower limit is constrained by accuracy of the tree code. To fulfil these requirements, we set Σ1 ≃ 0.1ΣO. 5.2 The dispersion relation To measure the dispersion relation, we compute Fourier transform of surface density for any frame of the simulation and then fit growth rate of each individual mode in time interval (⁠|$t^{\rm fit}_{\rm O}$|⁠, |$t^{{\rm fit}}_1)$|⁠. As the initial conditions do not correspond to the eigenfunctions, the modes relax at the beginning, and need some time before start growing at a temporarily constant growth rate. The choice of |$t^{\rm fit}_{\rm O}$| is constrained so as the modes with the longest wavelengths present in the computational domain already grow at a constant rate. |$t^{{\rm fit}}_1$| is constrained so as the modes with the highest growth rate have not reached the non-linear regime yet. We find that |$t^{\rm fit}_{\rm O}$| = 1/ωE78.FAST and |$t^{{\rm fit}}_1 = 3/\omega {_{\rm E78.FAST}}$| fulfil at best the requirements for our parameter choice. We calculate the dispersion relation for the same values of parameter A (models P18 and P99 in Table 2; right-hand panels of Fig. 2) as for the monochromatic models presented in Section 4. The data are binned to reduce noise. The critical wavenumber kMAX as well as the growth rate of wavenumbers k < kMAX is again in a very good agreement with E78 both in the self-gravity- and external pressure-dominated cases. Wavenumbers with k > kMAX are stable. We do not detect any significant deviation between our results and E78. 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, Jeans mass MJ, mean mass of gravitationally bound objects |$\bar{M}_{\rm CL}$| at fragmenting time tFRG, time tSINK, analytical e-folding time |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| for the most unstable wavenumber kFAST and transition time between linear and non-linear regime tNL. Run . A . nx × ny × nz . HHT/dz . . NJ . MJ . |$\bar{M}_{\rm CL}$| . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . tFRG . tSINK . tNL . . . . . . . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (Myr) . (Myr) . (Myr) . (Myr) . 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 Run . A . nx × ny × nz . HHT/dz . . NJ . MJ . |$\bar{M}_{\rm CL}$| . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . tFRG . tSINK . tNL . . . . . . . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (Myr) . (Myr) . (Myr) . (Myr) . 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 Open in new tab 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, Jeans mass MJ, mean mass of gravitationally bound objects |$\bar{M}_{\rm CL}$| at fragmenting time tFRG, time tSINK, analytical e-folding time |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| for the most unstable wavenumber kFAST and transition time between linear and non-linear regime tNL. Run . A . nx × ny × nz . HHT/dz . . NJ . MJ . |$\bar{M}_{\rm CL}$| . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . tFRG . tSINK . tNL . . . . . . . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (Myr) . (Myr) . (Myr) . (Myr) . 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 Run . A . nx × ny × nz . HHT/dz . . NJ . MJ . |$\bar{M}_{\rm CL}$| . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . tFRG . tSINK . tNL . . . . . . . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (Myr) . (Myr) . (Myr) . (Myr) . 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 Open in new tab 5.3 Evolution in non-linear regime Evolution of surface density for self-gravity-dominated model P99 is shown in Fig. 3. The fragmentation begins with emergence of round objects (at time from 0.7 to 2.1 Myr). These objects then gradually grow and become slender as was predicted by the second-order perturbation theory by Miyama et al. (1987b). The transformation of objects from roundish to filamentary-like is apparent between plots at 2.8 and 4.9 Myr. Sink particles form in the densest parts of the filaments and accrete material from surrounding filaments. Comparing the last plot at 4.9 Myr with a plot at the early stage of fragmentation (e.g. 2.1 Myr), we see that when a clump appears, it monolithically collapses to a gravitationally unstable object. Figure 3. Open in new tabDownload slide Evolution of surface density for the self-gravity-dominated layer (model P99). Number at the upper-left corner is the time in Myr. Sink particles are plotted as white circles. Area of a particular circle corresponds to the sink particle mass. Note that the colour scale for the upper and bottom rows is different. Figure 3. Open in new tabDownload slide Evolution of surface density for the self-gravity-dominated layer (model P99). Number at the upper-left corner is the time in Myr. Sink particles are plotted as white circles. Area of a particular circle corresponds to the sink particle mass. Note that the colour scale for the upper and bottom rows is different. Fragmentation in the pressure-dominated case is represented by model P18 (Fig. 4). At the beginning of fragmentation, the layer swiftly breaks into small objects (plots from 0.3 to 1.0 Myr) with masses smaller than the Jeans mass MJ. The sub-Jeans masses are a direct consequence of fragmentation according to E78 for layers with low A (right-hand panel of Fig. 1). For confinement of the objects, external pressure is more important than self-gravity, so apart from being immersed in an external gravitational field of the layer, the objects are equivalent to gravitationally stable Bonnor– Ebert spheres. The stable objects then gradually merge until enough mass for a gravitationally bound clump is assembled (see plots from 1.0 to 2.3 Myr). Merging often leads to non-radial accretion resulting in spinning up the fragments and disc formation around them (plot at time 2.3 Myr). As a bound clump is formed, it collapses and its cross-section for possible following mergers is reduced and merging rate decreases. Therefore, the fragmenting process in the pressure-dominated case, which proceeds via coalescence of many small clumps, is qualitatively different from the continuous collapse in the self-gravity-dominated case. We refer to the former and latter as coalescence-driven collapse and monolithic collapse, respectively. Figure 4. Open in new tabDownload slide Evolution of surface density for the external pressure-dominated layer (model P18). Caption is the same as in Fig. 3. Figure 4. Open in new tabDownload slide Evolution of surface density for the external pressure-dominated layer (model P18). Caption is the same as in Fig. 3. Since the analytical dispersion relations are often used for estimating fragmenting time-scales and mass of fragments, it is interesting to compare these quantities with those found in our simulations. We identify gravitationally bound fragments with the algorithm described in Appendix A. We experiment with several definitions of fragmenting time tFRG, and conclude that taking tFRG to be an instant when the total mass of gravitationally bound objects exceeds 1/2 of the total mass of the layer is a reasonable estimate for following reasons. As formation of bound objects starts, their total mass increases rapidly, so taking another fraction of the total mass than 1/2 would not lead to a significantly different time-scale. Time tFRG is also close to the time tSINK when sink particles in total contain 1/2 of the total mass of the layer (see Table 2 and left-hand panel of Fig. 5). Figure 5. Open in new tabDownload slide Fragmenting time and fragment mass of P models. Left-hand panel: time when 50 per cent of the total mass is in the form of sink particles (pluses) and gravitationally bound objects (circles). The solid line is the analytical estimate for E78, with e-folding time |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| scaled so as to match tFRG for model P99. Middle panel: number nEFOLD of e-folding times |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| needed to reach time tFRG. The solid line shows the fit (equation 19), and the dotted line shows the transition between linear and non-linear regime |$t{_{\rm NL}}/t_{{\rm E78.FAST}}^{{\rm EFOLD}}$|⁠. Right-hand panel: average mass of gravitationally bound objects at time tFRG (circles). Analytical estimates for E78, V83 and the Jeans mass for the mid-plane density are shown with solid, dotted and dashed lines, respectively. Figure 5. Open in new tabDownload slide Fragmenting time and fragment mass of P models. Left-hand panel: time when 50 per cent of the total mass is in the form of sink particles (pluses) and gravitationally bound objects (circles). The solid line is the analytical estimate for E78, with e-folding time |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| scaled so as to match tFRG for model P99. Middle panel: number nEFOLD of e-folding times |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| needed to reach time tFRG. The solid line shows the fit (equation 19), and the dotted line shows the transition between linear and non-linear regime |$t{_{\rm NL}}/t_{{\rm E78.FAST}}^{{\rm EFOLD}}$|⁠. Right-hand panel: average mass of gravitationally bound objects at time tFRG (circles). Analytical estimates for E78, V83 and the Jeans mass for the mid-plane density are shown with solid, dotted and dashed lines, respectively. Comparison between tFRG and E78 estimate |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}= 1/\omega {_{\rm E78.FAST}}$| (equation 11) as a function of parameter A is plotted in the left-hand panel of Fig. 5. We test the commonly adopted assumption that the fragmentation occurs at a constant number of analytic e-folding times |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| regardless of A, so we multiply |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| by a constant to match the data point for run P99. The assumption of constant number of e-folding times holds in the self-gravity-dominated case, but it underestimates the collapse time when external pressure dominates. The parameters describing the layer are ΣO, cS and PEXT. From dimensional analysis, it follows that any quantity with the dimension of time is given by \begin{equation} t = \frac{c{_{\rm S}}f(A)}{G \Sigma {_{\rm O}}}, \end{equation} (15) where f(A) is an unknown function of the parameter A. Since both |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| and tFRG are given by equation (15) (with presumably different description of f(A)), the number of e-folding times when fragmentation happens |$n{_{\rm EFOLD}}= t{_{\rm FRG}}/t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| is only a function of A for any layer. If we obtain nEFOLD(A) from our simulations, we can estimate tFRG for any layer because |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| is already known from equation (11). The dependence of nEFOLD on A is shown in the middle panel of Fig. 5 (solid line). Whereas nEFOLD is nearly constant in the self-gravity-dominated case, it strongly increases as the external pressure increases. Before drawing conclusions from this result, we should verify that it is not simply caused by the fact that self-gravity-dominated models already start with effectively higher initial perturbations. For any model, the transition to the non-linear regime at time tNL (when Σ1 = ΣO) occurs at almost constant number of |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| (Table 2 and middle panel of Fig. 5), so the initial perturbations are in this sense comparable among models with different A. Taking tNL instead of t = 0 as the time when the perturbing amplitudes are comparable would even emphasize the dependence of the number of |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| on A; when A is near unity, the fragmentation occurs soon after tNL, when A is small, it takes many e-folding times to reach tFRG. This behaviour reflects the two different fragmenting scenarios; when A is near unity, the clumps continue collapsing, 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 correctly because fragmentation does not occur at a constant number of e-folding times for various A. The mean mass of bound objects |$\bar{M}_{\rm CL}$| at the fragmenting time as a function of A is listed in Table 2 and shown in the right-hand panel of Fig. 5. Analytical estimates for V83 and E78 (MV83 and ME78, respectively) and the Jeans mass MJ for the mid-plane density [we use |$M{_{\rm J}}=4.26 \, A c_s^4/(G^2 \Sigma {_{\rm O}})$| from equations 13–33 in Spitzer (1978)] are plotted by lines. The fragment masses for E78 and V83 are estimated from the wavelength with the highest growth rate (i.e. M = πΣO(λFAST/2)2, where λFAST = 2π/kFAST). In the self-gravity-dominated case, the mass of fragments is in a good agreement with the E78 prediction. This result is a natural consequence of the monolithic collapse since the mass of fragments formed in the linear regime does not significantly change later on, during the non-linear collapse. The thin shell approximation systematically underestimates the fragment masses for A ≳ 0.6. In the external pressure-dominated case, the mass of fragments is of the order of the Jeans mass. It is a result of their formation process, since the fragments are initially sub-Jeans and stable until they assembled approximately the Jeans mass and then collapse. Therefore, the estimate based on E78 leads to too small fragment mass. In contrast, V83 overestimates fragment mass as it does not take into account the decrease of the Jeans mass with increasing external pressure. 6 LAYERS CONFINED BY THERMAL PRESSURE AND POSSIBLE PATTERN FORMATION IN SURFACE DENSITY Fuchs (1996) proposes that when fragmentation of an infinitely thin disc becomes non-linear, a triple of modes with wavevectors |$\Vert \boldsymbol {k}\Vert \simeq k{_{\rm FAST}}$| inclined at angles around |$60 ^\circ$| (i.e. in the Fourier space forming an equilateral triangle) has the highest growth rate due to interaction. Wünsch & Palouš (2001) semi-analytically find similar behaviour for the surface of a shell assuming V83. As a result, the modes create a hexagonal pattern in surface density. On the other hand, applying second-order perturbation theory, Miyama et al. (1987a) show that a fragmenting layer breaks into gradually slendering filaments with no signs of the hexagonal pattern. In this section, we test the Fuchs proposition by searching regular patterns in surface density in our polychromatic simulations. Since we find no evidence for any pattern, we perform three dedicated models to test whether pattern formation arises at all under very idealized circumstances. Name of the special models exploring possible pattern formation is in the form I|${\langle }A\rangle {}\_{}{\langle }N_w{\rangle }$|_[S|D], where the first two numbers after ‘I’ represent the parameter A and the number after the underscore the number Nw of amplified wavepackets. For models with suffix S, the wavepackets are amplified by a smooth function (wavepackets and the function are described in Section 6.1), while for models with suffix D, the wavepackets are degenerated to single modes. 6.1 Initial conditions for perturbations For the dedicated I models, the amplitudes |$\tilde{A}(k)$| of the perturbing wavevectors are first generated by the same method as for the polychromatic models (see Section 5.1 for description). Then, for models I99_3_S and I99_3_D, we choose three modes |$\boldsymbol {k}_i$| (i = 1, 2, 3) with |$\Vert \boldsymbol {k}_i\Vert \simeq k{_{\rm E78.FAST}}$| inclined at an angle around 60° and amplify amplitudes of all modes inside a circle at centre |$\boldsymbol {k}_i$| of radius dk by a bell-shaped curve. Thus, the perturbing amplitudes |$\tilde{A}_n$| are calculated as \begin{equation} \tilde{A}_n(k) = \left\lbrace \begin{array}{l@{\quad}l} A_{{\rm int}} \tilde{A}(k) \exp (-\Vert \boldsymbol {k}, \boldsymbol {k}_i \Vert ^2_2/\sigma _A^2), & \Vert \boldsymbol {k}, \boldsymbol {k}_i \Vert _2 \le \mathrm{d}k \\ \\ \tilde{A}(k), & \mathrm{otherwise} \\ \end{array}\right. \end{equation} (16) where |$\Vert \boldsymbol {k}, \boldsymbol {k}_i \Vert _2 = \sqrt{(k_x - k_{ix})^2 + (k_y - k_{iy})^2}$| and Aint and σA are parameters of the multiplying curve used in the particular model (Table 3). The only parameter differentiating models I99_3_D and I99_3_S is the radius dk. We experimented with amplifying amplitudes of only the selected wavevectors |$\boldsymbol {k}_i$| (model I99_3_D) and wavepackets inside non-zero radius dk centred on |$\boldsymbol {k}_i$| (model I99_3_S). The initial conditions for model I99_1_D are identical to that of model I99_3_D except that only one mode, |$\boldsymbol {k}_1$|⁠, is amplified. Table 3. Parameters for simulations designed to study mode interaction (models I). First four columns have the same meaning as the columns in Table 1. Parameters Aint, σA and dk characterize properties of the mode amplifying function, equation (16). Further, we list indices i of amplified modes |$\boldsymbol {k}_i$|⁠. Time |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| is as in Table 2. Run . A . nx × ny × nz . HHT/dz . Aint . σA . dk . i . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . . . . . . (pc−1) . (pc−1) . . (Myr) . I99_3_S 0.99 512 × 512 × 32 5.1 10.0 28 14 1,2,3 0.66 I99_3_D 0.99 512 × 512 × 32 5.1 10.0 28 0 1,2,3 0.66 I99_1_D 0.99 512 × 512 × 32 5.1 10.0 28 0 1 0.66 Run . A . nx × ny × nz . HHT/dz . Aint . σA . dk . i . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . . . . . . (pc−1) . (pc−1) . . (Myr) . I99_3_S 0.99 512 × 512 × 32 5.1 10.0 28 14 1,2,3 0.66 I99_3_D 0.99 512 × 512 × 32 5.1 10.0 28 0 1,2,3 0.66 I99_1_D 0.99 512 × 512 × 32 5.1 10.0 28 0 1 0.66 Open in new tab Table 3. Parameters for simulations designed to study mode interaction (models I). First four columns have the same meaning as the columns in Table 1. Parameters Aint, σA and dk characterize properties of the mode amplifying function, equation (16). Further, we list indices i of amplified modes |$\boldsymbol {k}_i$|⁠. Time |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| is as in Table 2. Run . A . nx × ny × nz . HHT/dz . Aint . σA . dk . i . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . . . . . . (pc−1) . (pc−1) . . (Myr) . I99_3_S 0.99 512 × 512 × 32 5.1 10.0 28 14 1,2,3 0.66 I99_3_D 0.99 512 × 512 × 32 5.1 10.0 28 0 1,2,3 0.66 I99_1_D 0.99 512 × 512 × 32 5.1 10.0 28 0 1 0.66 Run . A . nx × ny × nz . HHT/dz . Aint . σA . dk . i . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . . . . . . (pc−1) . (pc−1) . . (Myr) . I99_3_S 0.99 512 × 512 × 32 5.1 10.0 28 14 1,2,3 0.66 I99_3_D 0.99 512 × 512 × 32 5.1 10.0 28 0 1,2,3 0.66 I99_1_D 0.99 512 × 512 × 32 5.1 10.0 28 0 1 0.66 Open in new tab 6.2 Evolution of interacting modes Polychromatic models P18 and P99 (i.e. the ones with extreme values of parameter A in our sample) enter non-linear regime at time 0.42 and 2.90 Myr, respectively (Table 2). Inspecting corresponding plots in Figs 3 and 4 by eye, we do not see pattern formation at any stage of fragmenting process. We perform more rigorous analysis to find a firmer support for this statement. Let Fourier transform of the surface density be B(kx, ky). A regular pattern in the surface density would manifest itself as an anisotropy of function B. We look at the plane (kx, ky) as it is given in polar coordinates specified by radius |$\Vert \boldsymbol {k}\Vert$| and azimuthal angle α. To search for the anisotropy, we analyse function B evaluated on two different radially averaged subsets in the (kx, ky) plane. In the first case, B is evaluated on a thin annulus with radius |$\Vert \boldsymbol {k}\Vert = k{_{\rm E78.FAST}}$|⁠, i.e. |$B{_{\rm RAD.1}}(\alpha ) \equiv B(\Vert \boldsymbol {k}\Vert = k{_{\rm E78.FAST}}, \alpha )$|⁠. In the second case, the annulus is wider so it includes more modes around radius kE78.FAST, i.e. |$B{_{\rm RAD.2}}(\alpha ) \equiv B(k{_{\rm E78.FAST}}/2 \le \Vert \boldsymbol {k}\Vert \le 3k{_{\rm E78.FAST}}/2, \alpha )$|⁠. If two modes inclined at an angle |$\tilde{\alpha }$| have the highest amplitudes, function BRAD.J(α) (letter J stands for 1 or 2) has maxima separated by the angle |$\tilde{\alpha }$|⁠. Thus, the amplitudes of modes adjoined by the angle α are proportional to the value of azimuthal autocorrelation of BRAD.J(α), i.e. \begin{equation} c_J(\alpha ) = {{\int }_{0}^{\pi}} B_{\rm RAD.J}({\alpha} ^{\prime }) B_{\rm RAD.J}({\alpha + \alpha }^{\prime }) \,\mathrm{d}{\alpha} ^{\prime }. \end{equation} (17) Function cJ(α) detects interaction of modes inclined by any angle, not necessarily |$60 ^\circ$|⁠. When analysing the data, it is important to note that the autocorrelation function always attains its maximum at zero. The maxima due to two modes adjoined by an angle are local maxima of function cJ. We denote by |$c_{_{\rm {MAX.}J}}$| the highest local maximum of function cJ after its global maximum cJ(0). During the simulations, the value of cJ increases as amplitudes of modes increase. Since the interacting modes are predicted to grow faster than the rest of modes, we normalize |$c_{_{\rm {MAX.}J}}$| with the mean 〈cJ〉. Time dependence of |$c_{_{\rm {MAX.}1}}/\langle c{_{J}} \rangle$| and |$c_{_{\rm {MAX.}2}}/\langle c{_{2}} \rangle$| for models P18, P60 and P99 is shown in the upper panel of Fig. 6. Vertical bars indicate the transition to the non-linear regime. We detect no significant growth of |$c_{_{\rm {MAX.}J}}/\langle c{_{J}} \rangle$| after the models enter the non-linear regime indicating isotropy in the plane B(kx, ky) and thus no formation of a regular pattern. Figure 6. Open in new tabDownload slide Upper panel: evolution of azimuthal autocorrelation, equation (17), for models P18 (black line), P60 (blue line) and P99 (red line). Subscript characters 1 and 2 refer to two versions of function BRAD.: BRAD.1 and BRAD.2. Vertical bars denote transition from linear to non-linear regime for particular model as represented by its colour. Bottom panel: angular separation α between modes with |$c_{_{\rm {MAX.}J}}/\langle c{_{J}} \rangle$| at given time for models P18 (black symbols), P60 (blue symbols) and P99 (red symbols). The symbols are circles and asterisks in linear and non-linear regime of fragmentation for each model, respectively. Open and filled symbols represent analysis with BRAD.1 and BRAD.2, respectively. The angle 60|$^\circ$|⁠, which is expected to be adjoined by the triple of the highest growing modes in scenario proposed by Fuchs (1996), is denoted by the vertical solid line. Figure 6. Open in new tabDownload slide Upper panel: evolution of azimuthal autocorrelation, equation (17), for models P18 (black line), P60 (blue line) and P99 (red line). Subscript characters 1 and 2 refer to two versions of function BRAD.: BRAD.1 and BRAD.2. Vertical bars denote transition from linear to non-linear regime for particular model as represented by its colour. Bottom panel: angular separation α between modes with |$c_{_{\rm {MAX.}J}}/\langle c{_{J}} \rangle$| at given time for models P18 (black symbols), P60 (blue symbols) and P99 (red symbols). The symbols are circles and asterisks in linear and non-linear regime of fragmentation for each model, respectively. Open and filled symbols represent analysis with BRAD.1 and BRAD.2, respectively. The angle 60|$^\circ$|⁠, which is expected to be adjoined by the triple of the highest growing modes in scenario proposed by Fuchs (1996), is denoted by the vertical solid line. The bottom panel of Fig. 6 shows the angle α at which |$c_{_{\rm {MAX.}J}}/\langle c{_{J}} \rangle$| is attained. Data points representing linear and non-linear regime are distinguished by circles and asterisks, respectively. In the linear regime, there is no preferred angle adjoined by modes with the highest growth rate. This result is not surprising since the modes are predicted to grow independently of one another in the linear regime. If modes started interacting in the non-linear regime, asterisks would be clustered around a particular angle (⁠|$\alpha = 60 ^\circ$| is represented by a vertical line) with concomitant increase of |$c_{_{\rm {MAX.}J}}/\langle c{_{J}} \rangle$|⁠. Instead, the symbols show neither preferential clustering nor a rapid increase of |$c_{_{\rm {MAX.}J}}/\langle c{_{J}} \rangle$|⁠. The absence of both signposts strongly disfavours the possibility of regular pattern formation. Models I enable us to study the evolution of amplitudes ai of the individual modes |$\boldsymbol {k}_i$|⁠; this is plotted in Fig. 7. Models I99_3_S (black lines) and I99_3_D (blue lines) contain the triple of modes, which according to Fuchs (1996) should have an increased growth rate in the non-linear regime. Red lines show numerical solution to the second-order perturbation equations for the thin shell approximation (cf. equations 28 in Fuchs 1996) with the initial mode amplitudes of model I99_3_D.2 The vertical bars mark the transition from the linear to non-linear regime. Recall that the growth rate ω is the time derivative of the plotted functions. The growth rate obtained from our simulations does not follow the increased growth rate predicted by Fuchs. In contrast, the simulated growth rate even decreases in the non-linear regime. Figure 7. Open in new tabDownload slide Amplitude evolution of the three amplified modes for models I99_3_S (black line) and I99_3_D (blue line), and for the single amplified mode for model I99_1_D (magenta line). The mode of model I99_3_D that is identical to the amplified mode of model I99_1_D is shown by the dashed blue line. Numerical solution to the second-order perturbation equations proposed by Fuchs (1996) is plotted by red lines. The transition from the linear to non-linear regime is shown by vertical bars. Figure 7. Open in new tabDownload slide Amplitude evolution of the three amplified modes for models I99_3_S (black line) and I99_3_D (blue line), and for the single amplified mode for model I99_1_D (magenta line). The mode of model I99_3_D that is identical to the amplified mode of model I99_1_D is shown by the dashed blue line. Numerical solution to the second-order perturbation equations proposed by Fuchs (1996) is plotted by red lines. The transition from the linear to non-linear regime is shown by vertical bars. In the previous paragraph, we demonstrate that the possible interaction between the triple of modes does not increase their growth rate above ωE78.FAST. Do the triple of modes interact at least to some degree? To assess this question, we use the same initial conditions as for model I99_3_D, but we amplify only one mode, |$\boldsymbol {k}_1$| (model I99_1_D). Amplitude of this mode (magenta line in Fig. 7) evolves very close to that of the corresponding mode in model I99_3_D (blue dashed line). Moreover, the similar shape between these two curves suggests that the non-linear terms due to the two other modes are not important. Note that to study pattern formation, we investigate mode interaction only in the azimuthal direction α of the Fourier space. We do not investigate mode interaction in radial direction |$\Vert \boldsymbol {k}\Vert$|⁠, which apparently arises in the non-linear regime. Our findings do not contradict the mechanism proposed by Miyama et al. (1987a,b), which leads to randomly oriented filaments. 7 LAYERS ACCRETING FROM ONE SIDE, AND CONTAINED BY THERMAL PRESSURE ON THE OTHER In this section, we study layers accreting homogeneous medium from the upper surface and bounded by thermal pressure from the lower surface. We investigate their dispersion relation and subsequent evolution in the non-linear regime. These models approximate a part of a shell sweeping up the ambient medium on one surface and backed by thermal pressure on the other. As we show in Section 8 that shells around H ii regions fragment in the pressure-dominated case, we start the simulations in this case, setting A = 0.3 at the beginning of a simulation. The generic name of accreting models is in the form A|${\langle }A\rangle {}\_{}{\langle }\scr {M}{\rangle }{}\_{\langle }T{_{\rm AMB}}\rangle \_[N$||L]. First number after letter A represents the value of the parameter A at the beginning, followed by the Mach number |$\scr {M}$| of the accreting medium and the temperature TAMB of the ambient medium imposing the thermal pressure. The suffix L or N indicates whether the simulation terminates in linear or non-linear regime, respectively. For example, simulation A30_08_10_L treats a layer with A = 0.30, which accretes at velocity |$\scr {M} = 8$|⁠, and which is backed with an ambient medium of temperature 10 000 K. We use two kinds of models to address two different issues. Models with suffix L are used for studying evolution in the linear regime (the dispersion relation and flows inside the layer). For this purpose, it is sufficient to use computational domain in xy directions significantly (eight times) smaller than in corresponding non-accreting model P30. Consequently, we can afford to use two times higher resolution HHT/dz, even though a number of grid cells in x and y directions are four times smaller. Models with suffix N focus on evolution in the non-linear regime, where the priority is to include many Jeans masses inside the computational domain. For this purpose, we use the same size of the computational domain and the resolution as for model P30. 7.1 Initial conditions for perturbations Initial perturbations are generated by the same method as for the polychromatic models (Section 5.1). 7.2 Evolution in linear regime Numerically obtained dispersion relations in time interval |$(t^{{\rm fit}}_{{\rm O}}, t^{{\rm fit}}_1) = (1/\omega {_{\rm E78.FAST}}, 2/\omega {_{\rm E78.FAST}})$| for models with TAMB = 300 K, i.e. A30_08_03_L, A30_20_03_L and A30_50_03_L (Table 4), are plotted in the upper panel of Fig. 8. As the layer accretes, its surface density and parameter A increase, we list these quantities at |$t^{\rm fit}_{\rm O}$| and |$t^{{\rm fit}}_1$| in the last four columns of Table 4. We compare the results with E78 for two extreme values of parameter A attained in the simulations; one corresponds to t = 0, i.e. A = 0.30, the other to the model with the highest A at |$t^{{\rm fit}}_1$|⁠, i.e. A = 0.48. If the dispersion relation for the accreting layer were the same as for the thermal pressure-confined layer, all models would fit between these curves. However, the range of unstable wavenumbers extends to significantly higher values of kHHT. The highest growth rate is around a factor of 2 higher than that of E78. To study possible influence of the temperature in the ambient medium, we calculate the same models with temperature TAMB = 104 K (lower panel of Fig. 8). Both kinds of models reproduce similar features indicating that the dispersion relation depends only weakly on the temperature of the ambient medium. Note that these features are unlikely to be caused by resolution, which is here higher (7.8 cells per HHT; Table 4) than that for polychromatic models (≃4HHT; Table 2), which reproduce their analytic dispersion relation very well. Figure 8. Open in new tabDownload slide Dispersion relations for layers with A = 0.3 accreting homogeneous medium at Mach numbers 8, 20 and 50 (markers). Solid and dashed lines show E78 dispersion relation for two extreme values of A: A = 0.3 and A = 0.48, respectively. Upper panel: temperature of the ambient medium is 300 K (models A30_08_03_L, A30_20_03_L and A30_50_03_L). The dotted line is to guide the eye for |$\scr {M}$| = 8 model. Bottom panel: temperature of the ambient medium is 104 K (models A30_08_10_L, A30_20_10_L and A30_50_10_L). Figure 8. Open in new tabDownload slide Dispersion relations for layers with A = 0.3 accreting homogeneous medium at Mach numbers 8, 20 and 50 (markers). Solid and dashed lines show E78 dispersion relation for two extreme values of A: A = 0.3 and A = 0.48, respectively. Upper panel: temperature of the ambient medium is 300 K (models A30_08_03_L, A30_20_03_L and A30_50_03_L). The dotted line is to guide the eye for |$\scr {M}$| = 8 model. Bottom panel: temperature of the ambient medium is 104 K (models A30_08_10_L, A30_20_10_L and A30_50_10_L). Table 4. Parameters of accreting simulations intended to study evolution in the linear regime (models A). First four columns have the same meaning as the columns in Table 1 (parameter A is taken at time zero). Further we list Mach number |$\scr {M}$| and density ρACC of the accreting medium, temperature of the ambient medium TAMB, and surface density and the value of parameter A at time |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| and |$2t_{{\rm E78.FAST}}^{{\rm EFOLD}}$|⁠. Surface density at the beginning is denoted by ΣO. Run . A . nx × ny × nz . HHT/dz . |$\scr {M}$| . ρACC . TAMB . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . AT1 . AT2 . ΣT1/ΣO . ΣT2/ΣO . . . . . . (10−21 g cm−3) . (K) . . . . . . A30_08_10_L 0.30 128 × 128 × 64 7.8 8 1.660 10 000 0.25 0.39 0.47 1.37 1.76 A30_20_10_L 0.30 128 × 128 × 64 7.8 20 0.320 10 000 0.25 0.33 0.38 1.16 1.35 A30_50_10_L 0.30 128 × 128 × 64 7.8 50 0.052 10 000 0.25 0.31 0.33 1.06 1.15 A30_08_03_L 0.30 128 × 128 × 64 7.8 8 1.660 300 0.25 0.39 0.48 1.39 1.78 A30_20_03_L 0.30 128 × 128 × 64 7.8 20 0.320 300 0.25 0.34 0.41 1.17 1.48 A30_50_03_L 0.30 128 × 128 × 64 7.8 50 0.052 300 0.25 0.31 0.39 1.07 1.39 Run . A . nx × ny × nz . HHT/dz . |$\scr {M}$| . ρACC . TAMB . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . AT1 . AT2 . ΣT1/ΣO . ΣT2/ΣO . . . . . . (10−21 g cm−3) . (K) . . . . . . A30_08_10_L 0.30 128 × 128 × 64 7.8 8 1.660 10 000 0.25 0.39 0.47 1.37 1.76 A30_20_10_L 0.30 128 × 128 × 64 7.8 20 0.320 10 000 0.25 0.33 0.38 1.16 1.35 A30_50_10_L 0.30 128 × 128 × 64 7.8 50 0.052 10 000 0.25 0.31 0.33 1.06 1.15 A30_08_03_L 0.30 128 × 128 × 64 7.8 8 1.660 300 0.25 0.39 0.48 1.39 1.78 A30_20_03_L 0.30 128 × 128 × 64 7.8 20 0.320 300 0.25 0.34 0.41 1.17 1.48 A30_50_03_L 0.30 128 × 128 × 64 7.8 50 0.052 300 0.25 0.31 0.39 1.07 1.39 Open in new tab Table 4. Parameters of accreting simulations intended to study evolution in the linear regime (models A). First four columns have the same meaning as the columns in Table 1 (parameter A is taken at time zero). Further we list Mach number |$\scr {M}$| and density ρACC of the accreting medium, temperature of the ambient medium TAMB, and surface density and the value of parameter A at time |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| and |$2t_{{\rm E78.FAST}}^{{\rm EFOLD}}$|⁠. Surface density at the beginning is denoted by ΣO. Run . A . nx × ny × nz . HHT/dz . |$\scr {M}$| . ρACC . TAMB . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . AT1 . AT2 . ΣT1/ΣO . ΣT2/ΣO . . . . . . (10−21 g cm−3) . (K) . . . . . . A30_08_10_L 0.30 128 × 128 × 64 7.8 8 1.660 10 000 0.25 0.39 0.47 1.37 1.76 A30_20_10_L 0.30 128 × 128 × 64 7.8 20 0.320 10 000 0.25 0.33 0.38 1.16 1.35 A30_50_10_L 0.30 128 × 128 × 64 7.8 50 0.052 10 000 0.25 0.31 0.33 1.06 1.15 A30_08_03_L 0.30 128 × 128 × 64 7.8 8 1.660 300 0.25 0.39 0.48 1.39 1.78 A30_20_03_L 0.30 128 × 128 × 64 7.8 20 0.320 300 0.25 0.34 0.41 1.17 1.48 A30_50_03_L 0.30 128 × 128 × 64 7.8 50 0.052 300 0.25 0.31 0.39 1.07 1.39 Run . A . nx × ny × nz . HHT/dz . |$\scr {M}$| . ρACC . TAMB . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . AT1 . AT2 . ΣT1/ΣO . ΣT2/ΣO . . . . . . (10−21 g cm−3) . (K) . . . . . . A30_08_10_L 0.30 128 × 128 × 64 7.8 8 1.660 10 000 0.25 0.39 0.47 1.37 1.76 A30_20_10_L 0.30 128 × 128 × 64 7.8 20 0.320 10 000 0.25 0.33 0.38 1.16 1.35 A30_50_10_L 0.30 128 × 128 × 64 7.8 50 0.052 10 000 0.25 0.31 0.33 1.06 1.15 A30_08_03_L 0.30 128 × 128 × 64 7.8 8 1.660 300 0.25 0.39 0.48 1.39 1.78 A30_20_03_L 0.30 128 × 128 × 64 7.8 20 0.320 300 0.25 0.34 0.41 1.17 1.48 A30_50_03_L 0.30 128 × 128 × 64 7.8 50 0.052 300 0.25 0.31 0.39 1.07 1.39 Open in new tab Evolution of density and velocity field inside the accreting layer (model A30_08_10_L) is shown in Fig. 9. As the upper part of the layer is mildly corrugated, the different directions of incidence between the ram pressure and thermal pressure lead to a dynamical instability similar to the one described by Vishniac (1983), and it results in gas motions towards the convex parts (panel a). The gas then moves through the layer and protrudes from contact discontinuity until it is eventually stopped by gravitational attraction of the layer [panels (b) and (c)]. Formation of the protrusions at the contact discontinuity broadens the unstable range of the dispersion relation towards higher kHHT as is apparent from Fig. 8. At the same time, the corrugations of the bottom surface of the layer shield the hot gas in concave regions, leading to a pressure drop and thus cause flows of hot gas towards the concave regions. Upward flows of hot gas cause upward motions of cold gas inside the layer [panels (b) and (c)]. However, the different nature of BCs leads to different behaviour of the flows when they reach the surface. Whereas the upper flows are eroded by a shock, the downward flows can significantly corrugate the contact discontinuity (panel d) because this surface is pliable. The self-gravity then pulls the protruding gas, while it is still replenished from the original direction, resulting in circular motions (panel d). At this time (ωE78.FASTt ≃ 4), the growth of perturbations with kHHT ≳ 0.7 is saturated, so the range of unstable wavenumbers is comparable to that of the thermal pressure-confined layer. Figure 9. Open in new tabDownload slide Density profile and velocity field on cross-sections for model A30_08_10_L. The gas is accreted from the top. The velocity vectors are plotted only for the gas inside the layer to avoid confusion. The length of the thick white arrow in panel (a) corresponds to velocity 1 Mach number inside the layer. Figure 9. Open in new tabDownload slide Density profile and velocity field on cross-sections for model A30_08_10_L. The gas is accreted from the top. The velocity vectors are plotted only for the gas inside the layer to avoid confusion. The length of the thick white arrow in panel (a) corresponds to velocity 1 Mach number inside the layer. 7.3 Evolution in non-linear regime The subsequent non-linear fragmentation is investigated by models A30_05_03_N, A30_10_03_N and A30_20_03_N (Table 5). Evolution of surface density is shown in Fig. 10. Model A30_05_03_N accretes at the highest rate (cf. the sixth column in Table 5), and substantial part of fragmentation occurs in the self-gravity-dominated case (A ≳ 0.6; value of A is at the upper-right corner of each frame of Fig. 10). Fragments emerge at t = 1 Myr and subsequently collapse, forming filaments and then sink particles at their junctions (t = 1.5 and 2 Myr). Although the BCs at the upper surface are different from that of purely thermal pressure-confined layers, we see a similar evolutionary pattern as for the monolithic collapse (model P99 described in Section 5.3). Figure 10. Open in new tabDownload slide Evolution of surface density for accreting simulations. Top row: |$\scr {M}$| = 5 (model A30_05_03_N), middle: |$\scr {M}$| = 10 (A30_10_03_N), bottom: |$\scr {M}$| = 20 (A30_20_03_N). Time in Myr is in the upper-left corner, and the value of parameter A in the upper-right corner. Figure 10. Open in new tabDownload slide Evolution of surface density for accreting simulations. Top row: |$\scr {M}$| = 5 (model A30_05_03_N), middle: |$\scr {M}$| = 10 (A30_10_03_N), bottom: |$\scr {M}$| = 20 (A30_20_03_N). Time in Myr is in the upper-left corner, and the value of parameter A in the upper-right corner. Table 5. Parameters of accreting simulations intended for the non-linear regime of fragmentation (models A). First four columns have the same meaning as the columns in Table 1, and quantities NJ, MJ, |$\bar{M}_{\rm CL}$|⁠, tE78.FAST and tFRG are explained in Table 2. We further list the Mach number |$\scr {M}$| and density ρACC of the accreting medium, accreting rate |$\dot{\Sigma }$|⁠, and values of parameter A and surface density at the time tFRG (ΣO is the surface density at the beginning). Run . A . nx × ny × nz . |$\frac{H{_{\rm HT}}}{\mathrm{d}z}$| . |$\scr {M}$| . ρACC . |$\dot{\Sigma }$| . NJ . MJ . |$\bar{M}_{\rm CL}$| . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . tFRG . AFRG . |$\frac{\Sigma {_{\rm FRG}}}{\Sigma {_{\rm O}}}$| . . . . . . (10−21 g cm−3) . |$\frac{\, \mathrm{M}_{_{\odot }}}{\, \mathrm{Myr}\,\, \mathrm{pc}^2}$| . . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (Myr) . (Myr) . A30_05_03_N 0.30 512 × 512 × 64 3.9 5 4.94 76 9.8 3.6 4.0 0.25 1.5 0.73 3.5 A30_10_03_N 0.30 512 × 512 × 64 3.9 10 1.27 39 9.8 3.6 3.5 0.25 1.7 0.59 2.4 A30_20_03_N 0.30 512 × 512 × 64 3.9 20 0.32 20 9.8 3.6 3.4 0.25 1.8 0.59 2.4 Run . A . nx × ny × nz . |$\frac{H{_{\rm HT}}}{\mathrm{d}z}$| . |$\scr {M}$| . ρACC . |$\dot{\Sigma }$| . NJ . MJ . |$\bar{M}_{\rm CL}$| . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . tFRG . AFRG . |$\frac{\Sigma {_{\rm FRG}}}{\Sigma {_{\rm O}}}$| . . . . . . (10−21 g cm−3) . |$\frac{\, \mathrm{M}_{_{\odot }}}{\, \mathrm{Myr}\,\, \mathrm{pc}^2}$| . . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (Myr) . (Myr) . A30_05_03_N 0.30 512 × 512 × 64 3.9 5 4.94 76 9.8 3.6 4.0 0.25 1.5 0.73 3.5 A30_10_03_N 0.30 512 × 512 × 64 3.9 10 1.27 39 9.8 3.6 3.5 0.25 1.7 0.59 2.4 A30_20_03_N 0.30 512 × 512 × 64 3.9 20 0.32 20 9.8 3.6 3.4 0.25 1.8 0.59 2.4 Open in new tab Table 5. Parameters of accreting simulations intended for the non-linear regime of fragmentation (models A). First four columns have the same meaning as the columns in Table 1, and quantities NJ, MJ, |$\bar{M}_{\rm CL}$|⁠, tE78.FAST and tFRG are explained in Table 2. We further list the Mach number |$\scr {M}$| and density ρACC of the accreting medium, accreting rate |$\dot{\Sigma }$|⁠, and values of parameter A and surface density at the time tFRG (ΣO is the surface density at the beginning). Run . A . nx × ny × nz . |$\frac{H{_{\rm HT}}}{\mathrm{d}z}$| . |$\scr {M}$| . ρACC . |$\dot{\Sigma }$| . NJ . MJ . |$\bar{M}_{\rm CL}$| . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . tFRG . AFRG . |$\frac{\Sigma {_{\rm FRG}}}{\Sigma {_{\rm O}}}$| . . . . . . (10−21 g cm−3) . |$\frac{\, \mathrm{M}_{_{\odot }}}{\, \mathrm{Myr}\,\, \mathrm{pc}^2}$| . . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (Myr) . (Myr) . A30_05_03_N 0.30 512 × 512 × 64 3.9 5 4.94 76 9.8 3.6 4.0 0.25 1.5 0.73 3.5 A30_10_03_N 0.30 512 × 512 × 64 3.9 10 1.27 39 9.8 3.6 3.5 0.25 1.7 0.59 2.4 A30_20_03_N 0.30 512 × 512 × 64 3.9 20 0.32 20 9.8 3.6 3.4 0.25 1.8 0.59 2.4 Run . A . nx × ny × nz . |$\frac{H{_{\rm HT}}}{\mathrm{d}z}$| . |$\scr {M}$| . ρACC . |$\dot{\Sigma }$| . NJ . MJ . |$\bar{M}_{\rm CL}$| . |$t_{{\rm E78.FAST}}^{{\rm EFOLD}}$| . tFRG . AFRG . |$\frac{\Sigma {_{\rm FRG}}}{\Sigma {_{\rm O}}}$| . . . . . . (10−21 g cm−3) . |$\frac{\, \mathrm{M}_{_{\odot }}}{\, \mathrm{Myr}\,\, \mathrm{pc}^2}$| . . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (⁠|$\, \mathrm{M}_{_{\odot }}$|⁠) . (Myr) . (Myr) . A30_05_03_N 0.30 512 × 512 × 64 3.9 5 4.94 76 9.8 3.6 4.0 0.25 1.5 0.73 3.5 A30_10_03_N 0.30 512 × 512 × 64 3.9 10 1.27 39 9.8 3.6 3.5 0.25 1.7 0.59 2.4 A30_20_03_N 0.30 512 × 512 × 64 3.9 20 0.32 20 9.8 3.6 3.4 0.25 1.8 0.59 2.4 Open in new tab In contrast, the amount of gas accreted in model A30_20_03_N is substantially smaller and the main part of fragmentation is accomplished still in the pressure-dominated case. At the beginning of the fragmentation (t ≲ 1.0 Myr), the density inside the layer is approximately constant and the overdensities seen in panels at t = 0.5 and 1.0 Myr are corrugations of the contact discontinuity. The corrugations then gradually merge until they become gravitationally unstable and collapse (panel at t = 2.0 Myr). Note that the values of A seen in Fig. 10 are slightly higher than due to the accretion since cooling of hot gas at the contact discontinuity also contributes to the surface density, which increases A. The gradual merging forms less defined and more fluffy filaments than the direct collapse at smaller Mach numbers. Thus, the fragmenting process of model A30_20_03_N is very similar to the coalescence-driven collapse seen in model P18. The mean clump mass |$\bar{M}_{\rm CL}$| for models A30_05_03_N, A30_10_03_N and A30_20_03_N (3.4–|$4.0 \, \mathrm{M}_{_{\odot }}$|⁠, Table 5) is close to the Jeans mass (⁠|$3.6 \, \mathrm{M}_{_{\odot }}$| for A = 0.3). We interpret values of |$\bar{M}_{\rm CL}$| as follows. As the layer starts accreting in the pressure-dominated case, its density [according to equation (6)] and therefore the Jeans mass MJ is only a weak function of A. The highest change of MJ is for the model with the highest accreting rate, i.e. A30_05_03_N; MJ decreases only by a factor of 1.7 from beginning to the fragmenting time. Compared to non-accreting model P99, where |$\bar{M}_{\rm CL}= 2.1 \, \mathrm{M}_{_{\odot }}$| (Table 2), the accreting simulations have at most by a factor of 2 higher |$\bar{M}_{\rm CL}$|⁠. Therefore, the Jeans mass at t = 0 is a relatively good estimate for fragment mass also for accreting layers, so we can use it as a proxy when estimating fragments properties for an accreting shell in Section 8. Finally, we discuss the influence of our cooling method (described in Section 3.3) on the BCs. Recall that intermixing the warm gas into the layer and its imminent cooling deplete the warm gas on the lower side of the layer. The depleted gas is compensated by warm gas inflow, which causes ram pressure on the contact discontinuity. If the ram pressure is strong enough, it would change the contact discontinuity to a shock front that is an artificial and unwanted behaviour (e.g. surface-gravity waves would be suppressed), so we should check that this does not occur. From the beginning to the fragmenting time, the highest ratio between the ram and thermal pressure at the contact discontinuity is 0.03, 0.12 and 0.5 for models A30_05_03_N, A30_10_03_N and A30_20_03_N, respectively. Thus, the ram pressure is negligible when compared to the thermal pressure for models A30_05_03_N, A30_10_03_N and it is not dominant for model A30_20_03_N so that the contact discontinuity is preserved in all simulations. 8 FRAGMENTATION OF A SHELL SWEPT UP BY AN EXPANDING H ii REGION We use results of previous sections to estimate properties of a fragmenting shell (e.g. fragmenting time, mass) swept up around an H ii region. The H ii region is powered by a massive star emitting |$\dot{N}$| ionizing photons per second situated in a homogeneous medium with density ρACC (we use subscript ACC for the surrounding medium since it is accreted in the reference frame of the shell). For convenience, we define |$\dot{N}_{49} = [\dot{N} / 10^{49 }{\rm \,s^{-1}}]$|⁠, cs2 = [cS/0.2 km s−1] and n3 = [nACC/103 cm−3], where cS is sound velocity in the swept-up shell and nACC is the number density in the surrounding medium. Following Whitworth et al. (1994b), we assume that fragmentation is accomplished in the pressure-dominated case, which means that the shell fragments by coalescence-driven collapse. At the end, we show that our results are consistent with this assumption. 8.1 An analytic estimate for fragmentation of a shell swept up by an expanding H ii region The remarkable property of our accreting simulations is that their fragmenting time depends very weakly on the Mach number (Table 5) and differs at most by a factor of 1.5 from their non-accreting counterpart (model P30). Mass of the fragments also differs at most by a factor of 2 from the non-accreting model. This behaviour suggests that fragmenting time and mass of fragments of a non-accreting layer can be used to estimate corresponding quantities for an accreting layer. Further we assume that a thermal pressure-confined layer with time-dependent surface density Σ(t) and external pressure PEXT(t) fragments at any instant as the same layer with constant ΣO and PEXT. Thus, when fragmentation starts at time tO, it is accomplished at time t1 fulfilling \begin{equation} 1 = \int _{t{_{\rm O}}}^{t_1} \frac{\omega {_{\rm E78.FAST}}(t)}{n{_{\rm EFOLD}}(A(t))} \mathrm{d}t = \int _{t{_{\rm O}}}^{t_1} \frac{\sqrt{0.276} \pi G \Sigma (A(t))}{c{_{\rm S}}A(t) n{_{\rm EFOLD}}(A(t))} \mathrm{d}t, \end{equation} (18) where the second equality follows from equation (11). Fitting measured data (Table 2 and middle panel in Fig. 5), we obtain \begin{equation} n{_{\rm EFOLD}}= 1.63/A + 3.79. \end{equation} (19) We use our accreting models (with time-dependent Σ) to test the assumption made in the previous paragraph. Models A30_05_03_N, A30_10_03_N and A30_20_03_N fragment at time 1.5 (1.36) Myr, 1.7 (1.64) Myr and 1.8 (1.85) Myr, respectively, where the first number is obtained from the simulation and the bracketed number is from equation (18) showing a good agreement. The radius R of a shell around an H ii region expands according to R = RST(1 + 7cIIt/4RST)4/7 (Spitzer 1978), where RST and cII is Strömgren radius and sound velocity in the ionized gas, respectively. Surface density of the shell and ram pressure due to swept-up gas evolve as Σ = ρACCR(t)/3 and PRAM = ρACC(dR/dt)2. Thermal pressure PEXT inside the H ii region is equal to PRAM. Thus, approximating a small patch on the shell wall as a layer, its parameter A (equation 3) evolves as \begin{equation} A = \left\lbrace 1 + \frac{18 c_{{\rm II}}^2 }{\pi G \rho {_{\rm ACC}}(R{_{\rm ST}}(1 + 7c_{{\rm II}} t/(4R{_{\rm ST}})))^2 }\right\rbrace ^{-1/2}. \end{equation} (20) Evolution of PEXT, Σ and A for walls of shells powered by sources with |$\dot{N}_{49} = 0.1$| (approximately an O9V star; Martins, Schaerer & Hillier 2005), |$\dot{N}_{49} = 1$| (O6V star) and |$\dot{N}_{49} = 10$| (two O3V stars) expanding into a homogeneous ambient medium of particle density n3 = 0.1, n3 = 1 and n3 = 10, respectively, is shown in Fig. 11. Plotted models are extremes in the sense that any evolutionary path for a shell within a set |$\dot{N}_{49} \times n_{3}$| in range (0.1, 10) × (0.1, 10) lies in between them. The mentioned set |$\dot{N}_{49} \times n_{3}$| covers typical values for Galactic H ii regions. Time tO as determined from equation (22) below is plotted as the first asterisk on the evolutionary path. Following asterisks represent instants when integral in equation (18) reaches 0.25, 0.5, 0.75 and 1.0. They show that fragmentation proceeds in the pressure-dominated case. Figure 11. Open in new tabDownload slide Evolutionary paths for shell walls around expanding H ii regions. The paths with (⁠|$\dot{N}_{49}$|⁠, n3) = (0.1, 0.1), (⁠|$\dot{N}_{49}$|⁠, n3) = (1, 1) and (⁠|$\dot{N}_{49}$|⁠, n3) = (10, 10) are shown by red, black and blue lines, respectively. Circles with green numbers give time in Myr since the expansion started. White asterisks represent the instants when the integral in equation (18) equals 0, 0.25, 0.50, 0.75 and 1.0. Parameter A is shown by black labelled contours. Colour scale and dashed blue contours represent logarithm of the integrand in equation (18) expressed in Myr−1. The plot is constructed for cs2 = 1. Figure 11. Open in new tabDownload slide Evolutionary paths for shell walls around expanding H ii regions. The paths with (⁠|$\dot{N}_{49}$|⁠, n3) = (0.1, 0.1), (⁠|$\dot{N}_{49}$|⁠, n3) = (1, 1) and (⁠|$\dot{N}_{49}$|⁠, n3) = (10, 10) are shown by red, black and blue lines, respectively. Circles with green numbers give time in Myr since the expansion started. White asterisks represent the instants when the integral in equation (18) equals 0, 0.25, 0.50, 0.75 and 1.0. Parameter A is shown by black labelled contours. Colour scale and dashed blue contours represent logarithm of the integrand in equation (18) expressed in Myr−1. The plot is constructed for cs2 = 1. Substituting A from equation (20) for t in equation (18) gives \begin{equation} \frac{0.81 c{_{\rm S}}}{\left\lbrace (G \rho {_{\rm ACC}})^3 R{_{\rm ST}}^6 c_{{\rm II}}^8 \right\rbrace ^{1/14}} = \int _{A{_{\rm O}}}^{A_{1}} \frac{(1 - A^2)^{-25/14}}{n{_{\rm EFOLD}}(A) A^{3/7}} \mathrm{d}A, \end{equation} (21) where AO = A(tO) and A1 = A(t1) denote the value of parameter A at the beginning and end of fragmentation, respectively. We adopt tO as the earliest time when self-gravity of a fragment containing one Jeans mass (its radius is |$R{_{\rm J}}= \sqrt{M{_{\rm J}}/(\pi \Sigma )}$|⁠) dominates stretching, i.e. GΣ = RJ(4/(7tO))2 (see appendix A in Whitworth et al. 1994b).3 We neglect the pressure gradient term since the unstable objects form via coalescence when the role of pressure support is already lost. Thus, \begin{equation} t{_{\rm O}}= 0.76 \, \,{\rm Myr} \, c_{s2}^{\frac{28}{37}} \, n_3^{-\frac{33}{74}} \, \dot{N}_{49}^{-\frac{4}{37}}. \end{equation} (22) 8.2 The properties of fragments formed in a shell swept up by an expanding H ii region Integrating equation (21), we find A(t1) and from equation (20) fragmenting time t1. We evaluate integral equation (21) for cS = 0.2 km s−1 on a set of |$\dot{N}_{49} \times n_{3}$| in range (0.1, 10) × (0.1, 10), and by polynomial fitting we find \begin{eqnarray} t{_{\rm FRG}}&= 2.4 \, \, \mathrm{Myr}\, n_3^{-0.43} \, \dot{N}^{-0.12}_{49} \end{eqnarray} (23a) \begin{eqnarray} M{_{\rm FRG}}&= 3.1 \, \, \mathrm{M}_{_{\odot }}\, n_3^{-0.41} \, \dot{N}^{-0.16}_{49} \end{eqnarray} (23b) \begin{eqnarray} r{_{\rm FRG}}&= 0.12 \, \, \mathrm{pc}\, n_3^{-0.44} \, \dot{N}^{-0.11}_{49} \end{eqnarray} (23c) \begin{eqnarray} A{_{\rm FRG}}&= 0.51 \, n_3^{0.04} \, \dot{N}^{-0.08}_{49} \end{eqnarray} (23d) \begin{eqnarray} \Sigma {_{\rm FRG}}&= 63 \, \, \, \mathrm{M}_{_{\odot }}\, \, \mathrm{pc}^{-2}\, n_3^{0.46} \, \dot{N}^{0.08}_{49} \end{eqnarray} (23e) \begin{eqnarray} R{_{\rm FRG}}&= 7.7 \, {\rm pc} \, n_3^{-0.53} \, \dot{N}^{0.08}_{49}, \end{eqnarray} (23f) where MFRG is fragment mass, rFRG fragment radius, and AFRG, RFRG and ΣFRG are parameter A, shell radius and surface density at time tFRG, respectively. Since we show in Sections 5.3 and 7 that fragment mass of a pressure-confined layer is comparable to the Jeans mass, we use the Jeans mass as the estimate of fragment mass in equation (23b). We compare equation (23) with results of previous studies in Section 9.3. Throughout this section, we have assumed that the shell fragments in the external pressure-dominated case. The highest value attained by A in equation (23d) justifies this assumption (AFRG ≲ 0.6) unless n3 is either very high or |$\dot{N}_{49}$| very low. We restrict our study to sound velocity cs2 = 1 (i.e. shell temperature 10 K) to fulfil the condition AFRG ≲ 0.6 for realistic values of |$\dot{N}_{49}$| and n3 because the area in the |$(\dot{N}_{49}, n_3)$| plane where the condition is met decreases as cS increases. Since fragmentation of self-gravity-dominated shells is more complicated (time-dependent ρ0, gravitational focusing on emerging fragments, etc.), it is not clear what is the proxy for fragment mass in this case and we do not study it here. 9 DISCUSSION 9.1 Comparison of numerical with analytic dispersion relations The simulations in linear regime (when perturbed surface density Σ1 is smaller than unperturbed ΣO) clearly reproduce the semi-analytical estimate E78 in both the self-gravity- (A = 0.99) and pressure-dominated (A = 0.18) cases (Fig. 2). We do not detect any systematic difference from E78 confirming correctness of the derivation done by Elmegreen & Elmegreen (1978). Although W10 provides a good approximation to the range of unstable wavenumbers in the self-gravity-dominated case, it underestimates the highest unstable wavenumber kMAX for A ≲ 0.5 (Fig. 1). V83 overestimates kMAX by a factor of 2 for A = 1, predicts correct kMAX for A ≃ 0.6 and underestimates kMAX for lower A. Comparing to work of Dale et al. (2009) and Wünsch et al. (2010), our non-accreting models are more suitable for testing dispersion relations of a layer since they have time-independent ΣO and no stretching. They are also completely free of Vishniac and Rayleigh–Taylor instabilities. The lowest value of A that Wünsch et al. (2010) reach is 0.42 (bottom-right panel in their fig. 7). Although their results favour W10 to V83, they can hardly distinguish W10 from E78. We recover a less noisy dispersion relation and also reach lower value of A where we can detect the differences between W10 and E78. Since W10 was proposed to include the influence of the external pressure, it may be surprising why W10 differs from our results (and therefore E78) when A ≲ 0.5. We illustrate the limitation of W10 on a simple model where the layer is perturbed by an overdensity in a form of a circular patch of radius RC and central surface density Σ1. The perturbed surface density decreases to zero from the centre to the edge at RC. We estimate its stability against lateral collapse by comparing self-gravity with pressure gradient between the centre and the edge of the perturbation taken in the mid-plane. Then the condition for a lateral collapse is \begin{equation} G \beta \pi \Sigma _1 \gtrsim \frac{P_1 - P{_{\rm O}}}{\rho {_{\rm O}}R_{{\rm C}}}, \end{equation} (24) where P1 and PO is the pressure in the centre and the edge of the perturbation, respectively. Factor β takes into account that the perturbation generates non-spherically symmetric gravitational field. Factor β is of order unity. Substituting Σ1 + ΣO and ΣO into equation (6) and using equation (5), equation (24) becomes RCβ ≳ 2HHT. Assuming RC = λ/2, equation (24) takes the form \begin{equation} \beta \pi /2 \gtrsim k H{_{\rm HT}}. \end{equation} (25) Therefore, the model explains why the range of unstable wavelengths is of the order of the layer thickness for any A as proposed by E78 and found in our simulations. The limit of W10 (equation 9) in A → 0 predicts a significantly narrower range of unstable wavelengths, kHHT ≤ 9π2A2/20. The difference is due to how a perturbation is treated in the models. Imposing a perturbation with positive surface density Σ1 when A is small, the layer thickens locally, so the majority of perturbed mass is deposited above the surface with very little density enhancement in the mid-plane. The layer is almost of constant density where variations in surface density are dominated by corrugations of the surface. The weak dependence of P1 on Σ1 when 2PEXT/(πGΣO2) ≫ 1 is apparent from equation (6). The density (pressure) gradient is much smaller than it would be if the perturbation were placed in the mid-plane, so the pressure gradient is easily overcome by self-gravity and substantially shorter wavelengths are unstable. In contrast, this mechanism is absent in W10 where the perturbed mass is in the form of a homogeneous spheroid with no possibility to place the excess mass outside its surfaces. This assumption in W10 results in much higher difference between pressures PO and P1, which needs a more massive fragment (with longer λMAX) to be overcome with its self-gravity. In Section 8, assuming that shells swept up by H ii regions are at temperature 10 K, we conclude that the shells fragment in the pressure-dominated case. In this case, the significant difference between E78 and the other two dispersion relations leads to a substantially different prediction for fragment masses. On the other hand, V83 or W10 may still provide rough estimates for fragment masses of large H i shells, which are self-gravity dominated. Also note that E78 cannot be applied to layers forming at the interface between two colliding streams because such layers lack a surface with a contact discontinuity. 9.2 The effect of pre-existing density structure in the accretion flow We assume that the perturbing amplitudes for both polychromatic and accreting simulations do not depend on k (see Section 5.1 for details), and that the accreting simulations accrete homogeneous medium. These assumptions significantly constrain the possible spectrum of perturbations. This simple model was adopted in order to reduce the number of free parameters. The real ISM is highly inhomogeneous and when accreted, the resulting surface density enhancements have presumably different spectrum of perturbations. Since our fit giving the dependence of the number of e-folding times nEFOLD on A is done for the spectrum we adopt, the form of nEFOLD = nEFOLD(A) is presumably a function of the spectrum shape. In other words, the spectrum introduces another dimensionless parameters into equation (15). The parameters may further alter the properties of the fragmenting shell (equation 23). Accretion of even a single cloud could modify properties of final fragments substantially. For illustration, consider a spherically symmetric cloud. When accreted on a homogeneous layer, the ensuing structure would be a superposition of the layer and a circular overdensity. We further assume that the circular cloud contains many Jeans masses. This configuration is presumably prone to the edge effect (Burkert & Hartmann 2004; Pon et al. 2012), so it would tend to collapse laterally on its own time-scale tGC proportional to its radius rAC and overdensity Σ1 as |$t_{{\rm GC}} \simeq \sqrt{r_{{\rm AC}}/(\pi G \Sigma _1)}$|⁠. Given that the layer contains pre-existing perturbations of small amplitudes, the globally collapsing circular cloud would also fragment as a part of a layer and would break into pieces as we saw in Section 5.3. We suggest that external pressure could influence the result of fragmentation also in these circumstances if the circular cloud breaks into pieces before it collapsed globally. If the cloud is dominated by self-gravity, the fragments undergo an imminent collapse, so their radii shrink significantly diminishing probability of further coalescence. Consequently, the fragments would form a cluster in the centre. If the cloud is dominated by external pressure, the fragments are still stable against gravitational collapse, so their radii are comparable to their distances, and they can coalesce. However, the coalescence rate would be significantly enhanced in comparison to the scenario we saw in Section 5.3, since the fragments would be attracted to the centre or swept up by the passing edge of the cloud. This probably leads to objects substantially more massive than the Jeans mass. This issue may be an objective for future work. 9.3 Fragmenting shells in other works When estimating properties of fragments forming in the shell, we use significantly simplified models. The most important simplifications are the following: no deceleration, no ionizing radiation penetrating into the shell, no complicated thermodynamics and homogeneous medium into which the shell expands. Vertical acceleration of the unperturbed layer has reflective symmetry relative to the mid-plane; this assumption is violated in decelerating shells where it changes the dispersion relation as is shown by Iwasaki, Inutsuka & Tsuribe (2011a). The role of stretching is probably not crucial as estimated by Whitworth et al. (1994a). These phenomena may influence the course of fragmentation significantly and possibly modify or even suppress coalescence-driven collapse in real shells. Nevertheless, our work improves the model proposed by Whitworth et al. (1994a), which adopts the same simplified assumptions, since we use more appropriate estimates for mass of fragments and relevant time-scales. Taking into account approximations we made, we emphasize that equation (23) is rather an application of our results for layers than the definite answer to the problem of propagating star formation due to photoionizing feedback. Assuming temperature of 10 K in the shells, equation (23) predicts that shell fragmentation tends to form objects of too low mass for star formation to propagate. The dependence of our estimate equation (23) on |$\dot{N}_{49}$| and n3 is very close to the results provided by Whitworth et al. (1994a) and Iwasaki et al. (2011b, with sound velocity cs2 = 1). Compared to Whitworth et al. (1994a), equation (23a) predicts a slightly later fragmenting time ensuing in a larger shell radius (numerical constants they provide are 1.56 Myr and 5.8 pc). However, our model predicts much smaller fragment masses (3.1 versus 23 M⊙) and smaller fragment radii (0.12 versus 0.41 pc). This result is the consequence of taking as a proxy for fragment mass the Jeans mass for a significantly compressed layer instead of the Jeans mass for not compressed layer as done by Whitworth et al. (1994a). In order to test the numerical constants in formulae derived by Whitworth et al. (1994a), Dale, Bonnell & Whitworth (2007) select one point in the |$\dot{N}_{49}, n_3$| plane, |$(\dot{N}_{49}, n_3) = (1.0, 0.2)$| and follow its evolution by a smoothed particle hydrodynamics (SPH) simulation. They find an approximate agreement with the work of Whitworth et al. (1994a). However, comparing density of their shell ∼10−21 g cm−3 with the analytic value that is approximately 50 times higher, we suspect that their resolution is not enough to resolve the shell vertically. The low density in the shell leads to fragments of artificially high mass. On the other hand, our mass estimate in equation (23b) is close to the value 3.5 M⊙ found by Iwasaki et al. (2011b) based on monochromatic SPH simulations. We cannot compare fragmenting time since they use a different definition of the quantity. 10 CONCLUSIONS We perform 3D hydrodynamic simulations of self-gravitating isothermal layers in order to investigate their dispersion relations and subsequent fragmentation. We apply the results of fragmenting and accreting layers to estimate typical masses and fragmenting time for a shell swept up around an H ii region. We find that if the perturbations are small, thermal pressure-confined layers fragment in an excellent agreement with the dispersion relation proposed by Elmegreen & Elmegreen (1978) both in the self-gravity- and pressure-dominated cases. PAGI (Wünsch et al. 2010) is a good approximation and thin shell (Vishniac 1983) a rough approximation in the self-gravity-dominated case. However, the latter two are inappropriate for the pressure-dominated case (A ≲ 0.5; equation 3) because they do not take into account surface corrugations. Since the shells around H ii regions presumably fragment in the pressure-dominated case, we suggest that both PAGI and thin shell dispersion relations are not appropriate for their description. Fragmentation in the non-linear regime proceeds in two qualitatively different ways depending on parameter A. When the layer is dominated by self-gravity, fragments form by a monolithic collapse. When the layer is dominated by external pressure, it fragments in two steps; first it breaks into small gravitationally stable objects. It is an interesting feature because purely gravitational instability in planar geometry can form objects that are themselves stable against gravitational collapse. Then, the fragments continuously merge until they form an unstable object that eventually collapses (coalescence-driven collapse). Neither of the dispersion relations presented can predict correctly properties of gravitationally bound fragments for a pressure-dominated layer. However, a rough estimate for their mass is the Jeans mass in the mid-plane of the layer. In the non-linear regime, we investigate the possibility that the layer self-organizes and forms a regular pattern on its surface. We seek our standard models covering various degree of external pressure confinement (A = 0.18, A = 0.60, A = 0.99), and find no evidence for this scenario. In addition, we obtain the same result with two models designed to enhance pattern formation. We note that our method is able to find any regular pattern, not only a hexagonal pattern as is suggested to form by analytic work of Fuchs (1996). Instead of regular patterns, we observe formation of randomly oriented filamentary-like structures in agreement with Miyama et al. (1987a,b). For pressure-dominated layers, we substitute boundary conditions on one surface from thermal pressure for an accreting homogeneous medium. The dispersion relation for accreting layers has the range of unstable wavenumbers extended towards higher values of k in comparison to the dispersion relation for thermal pressure-confined layers. The highest growth rate is by a factor of ≃2 higher. In the non-linear regime of fragmentation, accreting layers also undergo monolithic or coalescence-driven collapse depending on the importance of confining pressure during the major part of fragmentation. Layers with high accreting rate become earlier dominated by self-gravity and collapse monolithically, while layers with lower accreting rate remain dominated by external pressure and fragment via coalescence-driven collapse. Fragmentation occurs at time comparable to thermal pressure-confined layers with the same instantaneous surface density and ambient pressure. The mass of gravitationally bound fragments is again comparable to the Jeans mass. We use our results from accreting layers to estimate fragment properties of a shell swept up around an expanding H ii region. For typical density (103 cm−3) and temperature (10 K) in molecular clouds, the fragmentation is accomplished while the shell is still dominated by the external pressure. This leads to fragment masses ≃ 3 M⊙. Stars formed of fragments of this mass are not able to ignite new H ii regions. This indicates that for star formation to propagate, either higher temperature in shells or a different scenario (e.g. geometry) is required. Acknowledgments The authors like to thank the referee, Sven Van Loo, for constructive comments that helped to significantly improve the paper. FD, RW and JP acknowledge support from the Czech Science Foundation grant P209/12/1795 and by the project RVO:67985815. The presented simulations were performed on IT4I facilities supported by the Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project ‘IT4Innovations National Supercomputing Center – LM2015070’. 1 " We distinguish a layer, which ideally is plane-parallel stratified, from a shell, which ideally is part, or the whole, of a hollow spherically symmetric structure. Although most of the observations motivating this work concern shells, the configurations we simulate in this paper are layers. Then we apply results derived for layers in the analysis of star formation triggered by expanding H ii regions – so in that analysis we also use the term shell, because the spatial curvature and velocity divergence of the shell can be neglected. 2 " Since the initial conditions for models I are not eigenfunctions, they do not exhibit the constant growth rate from the beginning. Equations (28) in Fuchs (1996) are derived for eigenfunctions, so they already starts with the constant growth rate. To compensate for the offset, we divide the latter by a constant so that their amplitudes are equal to that of the former at the end of the linear regime. 3 " Our choice of tO is an approximation since for smaller fragments with λMAX, λMAX < RJ self-gravity dominates stretching earlier than at tO. However, the results (equation 23) are weakly dependent on tO. Any numerical coefficients in equation (23) would differ by less than 20 per cent if we used tO = 0 instead of equation (22). REFERENCES Barnes J. , Hut P., 1986 , Nature , 324 , 446 Crossref Search ADS Bik A. et al. , 2010 , ApJ , 713 , 883 Crossref Search ADS Blaauw A. , 1964 , ARA&A , 2 , 213 Crossref Search ADS Blaauw A. , 1991 , Lada C. J., Kylafis N. D., The Physics of Star Formation and Early Stellar Evolution , Kluwer , Dordrecht , 125 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Boyd D. F. A. , Whitworth A. P., 2005 , A&A , 430 , 1059 Crossref Search ADS Brown A. G. A. , de Geus E. J., de Zeeuw P. T., 1994 , A&A , 289 , 101 Burkert A. , Hartmann L., 2004 , ApJ , 616 , 288 Crossref Search ADS Churchwell E. et al. , 2007 , ApJ , 670 , 428 Crossref Search ADS Colella P. , Woodward P. R., 1984 , J. Comput. Phys. , 54 , 174 Crossref Search ADS Dale J. E. , Bonnell I. A., Whitworth A. P., 2007 , MNRAS , 375 , 1291 Crossref Search ADS Dale J. E. , Wünsch R., Whitworth A., Palouš J., 2009 , MNRAS , 398 , 1537 Crossref Search ADS Dawson J. R. , Mizuno N., Onishi T., McClure-Griffiths N. M., Fukui Y., 2008 , MNRAS , 387 , 31 Crossref Search ADS Dawson J. R. , McClure-Griffiths N. M., Kawamura A., Mizuno N., Onishi T., Mizuno A., Fukui Y., 2011 , ApJ , 728 , 127 Crossref Search ADS Deharveng L. , Zavagno A., Caplan J., 2005 , A&A , 433 , 565 Crossref Search ADS Deharveng L. et al. , 2010 , A&A , 523 , A6 Crossref Search ADS Doroshkevich A. G. , 1980 , SvA , 24 , 152 Egorov O. V. , Lozinskaya T. A., Moiseev A. V., Smirnov-Pinchukov G. V., 2014 , MNRAS , 444 , 376 Crossref Search ADS Ehlerová S. , Palouš J., 2005 , A&A , 437 , 101 Crossref Search ADS Elmegreen B. G. , 1994 , ApJ , 427 , 384 Crossref Search ADS Elmegreen B. G. , 1998 , Woodward C. E., Shull J. M., Thronson H. A. Jr, ASP Conf. Ser. Vol. 148, Origins , Astron. Soc. Pac. , San Francisco , 150 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Elmegreen B. G. , Elmegreen D. M., 1978 , ApJ , 220 , 1051 , (E78) Crossref Search ADS Elmegreen B. G. , Lada C. J., 1977 , ApJ , 214 , 725 Crossref Search ADS Ewald P. P. , 1921 , Ann. Phys. , 369 , 253 Crossref Search ADS Federrath C. , Banerjee R., Clark P. C., Klessen R. S., 2010 , ApJ , 713 , 269 Crossref Search ADS Fryxell B. et al. , 2000 , ApJS , 131 , 273 Crossref Search ADS Fuchs B. , 1996 , MNRAS , 278 , 985 Crossref Search ADS Goldreich P. , Lynden-Bell D., 1965 , MNRAS , 130 , 97 Crossref Search ADS Iwasaki K. , Inutsuka S.-i., Tsuribe T., 2011a , ApJ , 733 , 16 Crossref Search ADS Iwasaki K. , Inutsuka S.-i., Tsuribe T., 2011b , ApJ , 733 , 17 Crossref Search ADS Kim J.-G. , Kim W.-T., Seo Y. M., Hong S. S., 2012 , ApJ , 761 , 131 Crossref Search ADS Klessen R. , 1997 , MNRAS , 292 , 11 Crossref Search ADS Lubow S. H. , Pringle J. E., 1993 , MNRAS , 263 , 701 Crossref Search ADS MacNeice P. , Olson K. M., Mobarry C., de Fainchtein R., Packer C., 2000 , Comput. Phys. Commun. , 126 , 330 Crossref Search ADS Martins F. , Schaerer D., Hillier D. J., 2005 , A&A , 436 , 1049 Crossref Search ADS Miyama S. M. , Narita S., Hayashi C., 1987a , Prog. Theor. Phys. , 78 , 1051 Crossref Search ADS Miyama S. M. , Narita S., Hayashi C., 1987b , Prog. Theor. Phys. , 78 , 1273 Crossref Search ADS Pon A. , Toalá J. A., Johnstone D., Vázquez-Semadeni E., Heitsch F., Gómez G. C., 2012 , ApJ , 756 , 145 Crossref Search ADS Salmon J. K. , Warren M. S., 1994 , J. Comput. Phys. , 111 , 136 Crossref Search ADS Simpson R. J. et al. , 2012 , MNRAS , 424 , 2442 Crossref Search ADS Spitzer L. , 1978 , Physical Processes in the Interstellar Medium , Wiley , New York Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Spitzer L. Jr, 1942 , ApJ , 95 , 329 Crossref Search ADS Springel V. , 2005 , MNRAS , 364 , 1105 Crossref Search ADS Truelove J. K. , Klein R. I., McKee C. F., Holliman J. H. II, Howell L. H., Greenough J. A., 1997 , ApJ , 489 , L179 Crossref Search ADS Van Loo S. , Keto E., Zhang Q., 2014 , ApJ , 789 , 37 Crossref Search ADS Vishniac E. T. , 1983 , ApJ , 274 , 152 , (V83) Crossref Search ADS Whitworth A. P. , Bhattal A. S., Chapman S. J., Disney M. J., Turner J. A., 1994a , A&A , 290 , 421 Whitworth A. P. , Bhattal A. S., Chapman S. J., Disney M. J., Turner J. A., 1994b , MNRAS , 268 , 291 Crossref Search ADS Wünsch R. , Palouš J., 2001 , A&A , 374 , 746 Crossref Search ADS Wünsch R. , Dale J. E., Palouš J., Whitworth A. P., 2010 , MNRAS , 407 , 1963 , (W10) Crossref Search ADS APPENDIX A: CLUMP FINDING ALGORITHM In order to determine the fragmenting time and mass of emerging bound objects, we need to identify gravitationally bound entities. We search for a gravitationally bound object around each minimum in gravitational potential in the computational domain. The gravitational potential is a sum of the gaseous potential as determined by flash code and the potential due to sink particles. When investigating a particular minimum i at potential ϕi (Fig. A1), we draw a sphere of the sink particle radius rSINK with centre at the potential minimum, and check whether the sphere contains at least one grid cell with negative binding energy eCELL. The binding energy of a cell is \begin{eqnarray} e{_{\rm CELL}}& = & \frac{1}{2} m{_{\rm CELL}}(v{_{\rm CELL}}- v_{{\rm CL}})^2 + \frac{3 m{_{\rm CELL}}cs{_{\rm CELL}}^2}{2} \nonumber \\ &&+\, \frac{1}{2} m{_{\rm CELL}}(\phi {_{\rm CELL}}- \phi ^{ij}), \end{eqnarray} (A1) where mCELL, ϕCELL, vCELL and csCELL are mass, potential, bulk velocity and sound velocity of the cell. The kinetic energy is corrected by velocity of the mass centre vCL of the bound object, which is composed of all bound grid cells. Figure A1. Open in new tabDownload slide Position of saddles and potential minima to illustrate our clump finding algorithm. Potential minima and saddle points are marked by subscripts and superscripts, respectively. Equipotential curves are plotted by thin solid lines, and cells identified as bound are plotted by colour corresponding to their density. Note that only selected minima and saddles are marked. Figure A1. Open in new tabDownload slide Position of saddles and potential minima to illustrate our clump finding algorithm. Potential minima and saddle points are marked by subscripts and superscripts, respectively. Equipotential curves are plotted by thin solid lines, and cells identified as bound are plotted by colour corresponding to their density. Note that only selected minima and saddles are marked. The gravitational binding energy is related to the potential at the saddle point ϕij (the second letter in the superscript denotes the deeper minimum, i.e. ϕj < ϕi). The saddle for minimum i is a pair of two neighbouring cells at the lowest potential ϕCELL where a path in potential from one cell monotonically descends to minimum i and a path from the other cell monotonically descends to another, deeper minimum j (see Fig. A1). Mixed BCs are taken into account, so the path can lead through a periodic boundary. No cell with ϕCELL > ϕij can be bound to the minimum i. If at least one bound cell is found, the radius of the sphere is increased by one grid cell size and condition (A1) is applied to any cell in the spherical shell between the new and previous searching radius. The procedure is repeated until the sphere exceeds the saddle point, or until it exceeds the cloud boundary, i.e. all cells in the spherical shell are unbound. If sink particles are present, they are very close to local potential minima and their mass is added to the mass budget of the bound object identified around the minimum. The computational domain typically contains many objects. Considerable number of the objects are interacting and overlapping. The algorithm first sorts the minima in ordering of decreasing potential and starts searching around the minimum with the highest value of potential and continues in direction of decreasing potential. This searching direction assures that the objects with minima at higher potential levels, which are usually situated close to more massive objects with minima at lower potential levels, are investigated before so that their mass is properly assigned to the lower mass objects. Any cell identified as bound is marked, and cannot be accessed later when investigating a deeper minimum. We illustrate some of the features of the algorithm in Fig. A1. Note that the ordering of potential minima (ϕm > ϕl > ϕi > ϕj > ϕk) is not reflected by ordering of saddles ϕlk > ϕjk > ϕmj > ϕij. Only cells around i with ϕ < ϕij can be assigned to minimum i. Since the saddle ϕjk for object j is at a higher equipotential surface than previous saddle ϕij, many cells situated close to i, which were not previously assigned to minimum i, are now identified as gravitationally bound and assigned to j. Cells situated around a shallow minimum at m, where no bound object was found previously, are checked if they are bound to j. The algorithm is able to find objects both filling and underfilling their critical equipotentials (e.g. i and l) and also identify highly non-spherically symmetric objects (e.g. j). © 2016 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
TI - Fragmentation of vertically stratified gaseous layers: monolithic or coalescence-driven collapse
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stw3354
DA - 2017-05-01
UR - https://www.deepdyve.com/lp/oxford-university-press/fragmentation-of-vertically-stratified-gaseous-layers-monolithic-or-67L2KBMkbN
SP - 4423
VL - 466
IS - 4
DP - DeepDyve
ER -