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

TL;DR

This paper uses 3D hydrodynamic simulations to study how vertically stratified gaseous layers fragment under different pressure conditions, revealing a dichotomy between monolithic collapse and coalescence-driven collapse.

Contribution

It identifies the key role of external pressure in determining the fragmentation mode of self-gravitating gaseous layers, providing new insights into star formation processes.

Findings

Layer fragmentation depends on external pressure dominance.

Coalescence occurs when external pressure dominates confinement.

Fragments are typically low-mass, limiting massive star formation.

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 (1978), and consequently deviates from the dispersion relations based on the thin shell approximation (Vishniac 1983) or pressure assisted gravitational instability (W\"unsch et al. 2010). 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 ( ~< 3 M_Sun ), and therefore they are unlikely to spawn stars that are sufficiently massive to promote sequential self-propagating star formation.

Tables

Table 1. Table 1: Parameters for simulations of layers confined by thermal pressure with monochromatic perturbations (models M). We list parameter A 𝐴 A , number of cells at the highest refinement level n x subscript 𝑛 𝑥 n_{x} , n y subscript 𝑛 𝑦 n_{y} and n z subscript 𝑛 𝑧 n_{z} , resolution in the vertical direction H HT / d ​ z subscript 𝐻 HT 𝑑 𝑧 H_{{}_{\rm HT}}/dz , parameter k ​ H HT 𝑘 subscript 𝐻 HT kH_{{}_{\rm HT}} where k 𝑘 k is the selected wavenumber, t E78 EFOLD superscript subscript 𝑡 E78 EFOLD t_{{}_{\rm E78}}^{{}_{\rm EFOLD}} is the analytical e–folding time according to E78 for given k 𝑘 k , and t NUM subscript 𝑡 NUM t_{{}_{\rm NUM}} ( t NUM = 1 / ω subscript 𝑡 NUM 1 𝜔 t_{{}_{\rm NUM}}=1/\omega ) is the e–folding time measured in our simulations.

Run	$A$	$n_x\times n_y\times n_z$	$H_{{}_{\mathrm{HT}}}/dz$	$k{H}_{{}_{\mathrm{HT}}}$	$t_{{}_{\mathrm{E78}}}^{{}_{\mathrm{EFOLD}}}$	$t_{{}_{\mathrm{NUM}}}$
					[Myr]	[Myr]
M18_020	0.18	640 $\times$ 128 $\times$ 128	20.2	0.20	0.16	0.16
M18_025	0.18	512 $\times$ 128 $\times$ 128	20.2	0.25	0.15	0.16
M18_033	0.18	384 $\times$ 128 $\times$ 128	20.2	0.33	0.15	0.16
M18_050	0.18	256 $\times$ 128 $\times$ 128	20.2	0.50	0.19	0.19
M18_062	0.18	208 $\times$ 128 $\times$ 128	20.2	0.62	0.41	0.41
M18_073	0.18	176 $\times$ 128 $\times$ 128	20.2	0.73	-	-
M18_089	0.18	144 $\times$ 128 $\times$ 128	20.2	0.89	-	-
M99_025	0.99	512 $\times$ 128 $\times$ 128	20.2	0.25	0.74	0.76
M99_033	0.99	384 $\times$ 128 $\times$ 128	20.2	0.33	0.69	0.70
M99_050	0.99	256 $\times$ 128 $\times$ 128	20.2	0.50	0.67	0.68
M99_073	0.99	176 $\times$ 128 $\times$ 128	20.2	0.73	0.77	0.78
M99_089	0.99	144 $\times$ 128 $\times$ 128	20.2	0.89	1.12	1.14
M99_114	0.99	112 $\times$ 128 $\times$ 128	20.2	1.14	-	-
M99_133	0.99	96 $\times$ 128 $\times$ 128	20.2	1.33	-	-

Table 2. 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 N J subscript 𝑁 J N_{{}_{\rm J}} , Jeans mass M J subscript 𝑀 J M_{{}_{\rm J}} , mean mass of gravitationally bound objects M ¯ CL subscript ¯ 𝑀 CL \bar{M}_{\rm CL} at fragmenting time t FRG subscript 𝑡 FRG t_{{}_{\rm FRG}} , time t SINK subscript 𝑡 SINK t_{{}_{\rm SINK}} , analytical e–folding time t E78 . FAST EFOLD superscript subscript 𝑡 formulae-sequence E78 FAST EFOLD t_{{}_{\rm E78.FAST}}^{{}_{\rm EFOLD}} for the most unstable wavenumber k FAST subscript 𝑘 FAST k_{{}_{\rm FAST}} and transition time between linear and non–linear regime t NL subscript 𝑡 NL t_{{}_{\rm NL}} .

Run	$A$	$n_x\times n_y\times n_z$	$H_{{}_{\mathrm{HT}}}/dz$		$N_{{}_J}$	$M_{{}_J}$	${\overline{M}}_{\mathrm{CL}}$	$t_{{}_{\mathrm{E78}.\mathrm{FAST}}}^{{}_{\mathrm{EFOLD}}}$	$t_{{}_{\mathrm{FRG}}}$	$t_{{}_{\mathrm{SINK}}}$	$t_{{}_{\mathrm{NL}}}$
						[$M_{{}_{\odot }}$]	[$M_{{}_{\odot }}$]	[Myr]	[Myr]	[Myr]	[Myr]
P18	0.18	1024 $\times$ 1024 $\times$ 64	4.0		8.8	2.3	1.1	0.15	2.01	2.15	0.42
P20	0.20	1024 $\times$ 1024 $\times$ 64	4.4		9.3	2.5	1.9	0.17	2.19	2.40	0.50
P22	0.22	1024 $\times$ 1024 $\times$ 64	4.4		12.8	2.8	1.7	0.18	2.11	2.23	0.58
P25	0.25	1024 $\times$ 1024 $\times$ 64	4.4		19.0	3.1	3.1	0.21	1.90	2.20	0.80
P30	0.30	512 $\times$ 512 $\times$ 64	3.9		9.8	3.6	2.1	0.25	2.20	2.40	0.99
P40	0.40	512 $\times$ 512 $\times$ 64	4.7		17.0	5.0	3.5	0.33	2.38	2.64	1.45
P60	0.60	512 $\times$ 512 $\times$ 64	6.6		29.7	7.4	10.0	0.47	3.00	3.35	2.22
P80	0.80	512 $\times$ 512 $\times$ 64	5.9		89.0	9.9	20.6	0.59	3.60	4.00	2.31
P99	0.99	512 $\times$ 512 $\times$ 64	9.1		71.0	12.5	32.0	0.66	4.15	4.90	2.90

Table 3. 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 A i ​ n ​ t subscript 𝐴 𝑖 𝑛 𝑡 A_{int} , σ A subscript 𝜎 𝐴 \sigma_{A} and d ​ k 𝑑 𝑘 dk characterise properties of the mode amplifying function, Eq. ( 16 ). Further, we list indices i 𝑖 i of amplified modes k → i subscript → 𝑘 𝑖 \vec{k}_{i} . Time t E78 . FAST EFOLD superscript subscript 𝑡 formulae-sequence E78 FAST EFOLD t_{{}_{\rm E78.FAST}}^{{}_{\rm EFOLD}} is as in Table 2 .

Run	$A$	$n_x\times n_y\times n_z$	$H_{{}_{\mathrm{HT}}}/dz$	$A_{int}$	${\sigma }_A$	$dk$	$i$	$t_{{}_{\mathrm{E78}.\mathrm{FAST}}}^{{}_{\mathrm{EFOLD}}}$
					$\left[{\mathrm{pc}}^{-1}\right]$	$\left[{\mathrm{pc}}^{-1}\right]$		[Myr]
I99_3_S	0.99	512 $\times$ 512 $\times$ 32	5.1	10.0	28	14	1,2,3	0.66
I99_3_D	0.99	512 $\times$ 512 $\times$ 32	5.1	10.0	28	0	1,2,3	0.66
I99_1_D	0.99	512 $\times$ 512 $\times$ 32	5.1	10.0	28	0	1	0.66

Table 4. 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 𝐴 A is taken at time zero). Further we list Mach number ℳ ℳ \mathscr{M} and density ρ ACC subscript 𝜌 ACC \rho_{{}_{\rm ACC}} of the accreting medium, temperature of the ambient medium T AMB subscript 𝑇 AMB T_{{}_{\rm AMB}} , and surface density and the value of parameter A 𝐴 A at time t E78 . FAST EFOLD superscript subscript 𝑡 formulae-sequence E78 FAST EFOLD t_{{}_{\rm E78.FAST}}^{{}_{\rm EFOLD}} and 2 ​ t E78 . FAST EFOLD 2 superscript subscript 𝑡 formulae-sequence E78 FAST EFOLD 2t_{{}_{\rm E78.FAST}}^{{}_{\rm EFOLD}} . Surface density at the beginning is denoted Σ O subscript Σ O \Sigma_{{}_{\rm O}} .

Run	$A$	$n_x\times n_y\times n_z$	$H_{{}_{\mathrm{HT}}}/dz$	$ℳ$	${\rho }_{{}_{\mathrm{ACC}}}$	$T_{{}_{\mathrm{AMB}}}$	$t_{{}_{\mathrm{E78}.\mathrm{FAST}}}^{{}_{\mathrm{EFOLD}}}$	$A_{T1}$	$A_{T2}$	${\Sigma }_{T1}/{\Sigma }_{{}_O}$	${\Sigma }_{T2}/{\Sigma }_{{}_O}$
					$\left[{10}^{-21}g.{\mathrm{cm}}^{-3}\right]$	[$K$]
A30_08_10_L	0.30	128 $\times$ 128 $\times$ 64	7.8	8	1.660	10000	0.25	0.39	0.47	1.37	1.76
A30_20_10_L	0.30	128 $\times$ 128 $\times$ 64	7.8	20	0.320	10000	0.25	0.33	0.38	1.16	1.35
A30_50_10_L	0.30	128 $\times$ 128 $\times$ 64	7.8	50	0.052	10000	0.25	0.31	0.33	1.06	1.15
A30_08_03_L	0.30	128 $\times$ 128 $\times$ 64	7.8	8	1.660	300	0.25	0.39	0.48	1.39	1.78
A30_20_03_L	0.30	128 $\times$ 128 $\times$ 64	7.8	20	0.320	300	0.25	0.34	0.41	1.17	1.48
A30_50_03_L	0.30	128 $\times$ 128 $\times$ 64	7.8	50	0.052	300	0.25	0.31	0.39	1.07	1.39

Table 5. 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 N J subscript 𝑁 J N_{{}_{\rm J}} , M J subscript 𝑀 J M_{{}_{\rm J}} , M ¯ CL subscript ¯ 𝑀 CL \bar{M}_{\rm CL} , t E78 . FAST subscript 𝑡 formulae-sequence E78 FAST t_{{}_{\rm E78.FAST}} and t FRG subscript 𝑡 FRG t_{{}_{\rm FRG}} are explained in Table 2 . We further list the Mach number ℳ ℳ \mathscr{M} and density ρ ACC subscript 𝜌 ACC \rho_{{}_{\rm ACC}} of the accreting medium, accreting rate Σ ˙ ˙ Σ \dot{\Sigma} , and values of parameter A 𝐴 A and surface density at the time t FRG subscript 𝑡 FRG t_{{}_{\rm FRG}} ( Σ O subscript Σ O \Sigma_{{}_{\rm O}} is the surface density at the beginning).

Run	$A$	$n_x\times n_y\times n_z$	$\frac{H_{{}_{\mathrm{HT}}}}{dz}$	$ℳ$	${\rho }_{{}_{\mathrm{ACC}}}$	$\dot{\Sigma }$	$N_{{}_J}$	$M_{{}_J}$	${\overline{M}}_{\mathrm{CL}}$	$t_{{}_{\mathrm{E78}.\mathrm{FAST}}}^{{}_{\mathrm{EFOLD}}}$	$t_{{}_{\mathrm{FRG}}}$	$A_{{}_{\mathrm{FRG}}}$	$\frac{{\Sigma }_{{}_{\mathrm{FRG}}}}{{\Sigma }_{{}_O}}$
					[${10}^{-21}g.{\mathrm{cm}}^{-3}$]	$\frac{M_{{}_{\odot }}}{\mathrm{Myr}{\mathrm{pc}}^2}$		[$M_{{}_{\odot }}$]	[$M_{{}_{\odot }}$]	[Myr]	[Myr]
A30_05_03_N	0.30	512 $\times$ 512 $\times$ 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 $\times$ 512 $\times$ 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 $\times$ 512 $\times$ 64	3.9	20	0.32	20	9.8	3.6	3.4	0.25	1.8	0.59	2.4

Equations

$$\rho (z)=\frac{{\rho }_{{}_{\mathrm{O}}}}{{\cosh }^2(z/H_{{}_{\mathrm{O}}})},$$

$$H_{{}_{\mathrm{O}}}=\frac{c_{{}_{\mathrm{S}}}}{\sqrt{2\pi G{\rho }_{{}_{\mathrm{O}}}}},$$

$$A=\frac{1}{\sqrt{1+\left(2P_{{}_{\mathrm{EXT}}}/\pi G{\Sigma }_{{}_{\mathrm{O}}}^2\right)}},$$

$$z_{{}_{\mathrm{MAX}}}=H_{{}_{\mathrm{O}}}{\mathrm{arcosh}}^{-1}\left(\frac{{\rho }_{{}_{\mathrm{O}}}{c}_{{}_{\mathrm{S}}}^2}{P_{{}_{\mathrm{EXT}}}}\right)=H_{{}_{\mathrm{O}}}{\mathrm{arcosh}}^{-1}\left(\frac{1}{\sqrt{1-A^2}}\right).$$

$$H_{{}_{\mathrm{HT}}}=\frac{{\Sigma }_{{}_{\mathrm{O}}}{c}_{{}_{\mathrm{S}}}^2}{2P_{{}_{\mathrm{EXT}}}+\pi G{\Sigma }_{{}_{\mathrm{O}}}^2}=\frac{c_{{}_{\mathrm{S}}}^2{A}^2}{\pi G{\Sigma }_{{}_{\mathrm{O}}}},$$

$${\rho }_{{}_{\mathrm{O}}}=\frac{2P_{{}_{\mathrm{EXT}}}+\pi G{\Sigma }_{{}_{\mathrm{O}}}^2}{2c_{{}_{\mathrm{S}}}^2}=\frac{\pi G{\Sigma }_{{}_{\mathrm{O}}}^2}{2c_{{}_{\mathrm{S}}}^2{A}^2}.$$

$$\text{mn@boldsymbol}{\omega }_{{}_{V83}}^2(k)=2\pi G{\Sigma }_{{}_{\mathrm{O}}}k-c_{{}_{\mathrm{S}}}^2{k}^2.$$

$$\ddot{r}$$

$$\text{mn@boldsymbol}{\omega }_{{}_{W10}}^2(k)$$

$$\text{mn@boldsymbol}{\omega }_{{}_{E78}}^2=2\pi G{\Sigma }_{{}_{\mathrm{O}}}k(\eta F_{{}_{\mathrm{J}}}+(1-\eta )F_{{}_{\mathrm{D}}})-c_{{}_{\mathrm{EFF}}}^2{k}^2.$$

$${\omega }_{{}_{E78.FAST}}^2=0.276\times 2\pi G{\rho }_{{}_{\mathrm{O}}}=0.276{\left(\frac{\pi G{\Sigma }_{{}_{\mathrm{O}}}}{c_{{}_{\mathrm{S}}}A}\right)}^2.$$

$$\Delta a_{(\mathrm{p})}^{\max }=\frac{GM}{d^2}{\left(\frac{h}{d}\right)}^{p+1},$$

$${\rho }_{{}_{\mathrm{AMB}}}(z)$$

$$\sigma (z)=\underset{z_{{}_{\mathrm{BOT}}}}{\overset{z}{\int }}\rho (z^{\prime })d{z}^{\prime }.$$

$$t=\frac{c_{{}_{\mathrm{S}}}f(A)}{G{\Sigma }_{{}_{\mathrm{O}}}},$$

$${\tilde{A}}_n(k)=\{\begin{array}{ll}{A}_{int}\tilde{A}(k)\exp (-\Vert \vec{k},{\vec{k}}_i{\Vert }_2^2/{\sigma }_A^2)& ,\Vert \vec{k},{\vec{k}}_i{\Vert }_2\le \mathrm{d}k\\ \tilde{A}(k)& ,\mathrm{otherwise}\end{array}$$

$$c_J(\alpha )={\int }_0^{\pi }{B}_{{}_{RAD.J}}({\alpha }^{\prime })B_{{}_{RAD.J}}(\alpha +{\alpha }^{\prime })\mathrm{d}{\alpha }^{\prime }.$$

$$1={\int }_{t_{{}_{\mathrm{O}}}}^{t_1}\frac{{\omega }_{{}_{E78.FAST}}(t)}{n_{{}_{\mathrm{EFOLD}}}(A(t))}\mathrm{d}t={\int }_{t_{{}_{\mathrm{O}}}}^{t_1}\frac{\sqrt{0.276}\pi G\Sigma (A(t))}{c_{{}_{\mathrm{S}}}A(t)n_{{}_{\mathrm{EFOLD}}}(A(t))}\mathrm{d}t,$$

$$n_{{}_{\mathrm{EFOLD}}}=1.63/A+\mathrm{3.79.}$$

$$A={\left\{1+\frac{18c_{II}^2}{\pi G{\rho }_{{}_{\mathrm{ACC}}}(R_{{}_{\mathrm{ST}}}(1+7c_{II}t/(4R_{{}_{\mathrm{ST}}})){)}^2}\right\}}^{-1/2}.$$

$$\frac{0.81c_{{}_{\mathrm{S}}}}{{\left\{(G{\rho }_{{}_{\mathrm{ACC}}}{)}^3{R}_{{}_{\mathrm{ST}}}^6{c}_{II}^8\right\}}^{1/14}}={\int }_{A_{{}_{\mathrm{O}}}}^{A_1}\frac{(1-A^2{)}^{-25/14}}{n_{{}_{\mathrm{EFOLD}}}(A)A^{3/7}}\mathrm{d}A,$$

$$t_{{}_{\mathrm{O}}}=0.76\mathrm{Myr}{c}_{s2}^{\frac{28}{37}}{n}_3^{-\frac{33}{74}}{\dot{N}}_{49}^{-\frac{4}{37}}.$$

$$t_{{}_{\mathrm{FRG}}}$$

$$M_{{}_{\mathrm{FRG}}}$$

$$r_{{}_{\mathrm{FRG}}}$$

$$A_{{}_{\mathrm{FRG}}}$$

$${\Sigma }_{{}_{\mathrm{FRG}}}$$

$$R_{{}_{\mathrm{FRG}}}$$

$$G\beta \pi {\Sigma }_1\text{ga}\frac{P_1-P_{{}_{\mathrm{O}}}}{{\rho }_{{}_{\mathrm{O}}}{R}_{{}_C}},$$

$$\beta \pi /2\text{ga}k{H}_{{}_{\mathrm{HT}}}.$$

$$e_{{}_{\mathrm{CELL}}}$$

Full text

††pagerange: Fragmentation of vertically stratified gaseous layers: monolithic or coalescence–driven collapse–A††pubyear: 2017

1 Introduction

Massive stars ($M_{\star}\ga 8\,\mathrm{M}_{{}_{\odot}}$) strongly influence the interstellar medium (ISM) surrounding them, mainly via photoionisation, stellar winds and supernova explosions. Elmegreen & Lada (1977) propose a mechanism (Collect and Collapse) whereby an over-pressured 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 ionisation front. The resulting shell of cool 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% of the infrared shells identified by Churchwell et al. (2007) encircle H ii regions ionised by O- and early B-type stars, suggesting that the shells are due to feedback from these stars. After a careful examination, Deharveng et al. (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 et al., 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 111We 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. , 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 et al. (1987a, b) extend the linearised 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 which 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 et al. (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 et al. (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 onto 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 organised 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 summarises our main conclusions. Appendix A describes our algorithm developed for finding gravitationally bound fragments.

2 Linearised theory of layer fragmentation

2.1 The unperturbed state

Before introducing the dispersion relations describing the growth rate of perturbations in the self-gravitating presure-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 $c_{{}_{\rm S}}$, so the density distribution is as derived by Spitzer (1942) and Goldreich & Lynden-Bell (1965), i.e.

[TABLE]

where $\rho_{{}_{\rm O}}$ is the density on the midplane ($z\!=\!0$),

[TABLE]

is the vertical scale height, and $G$ is the gravitational constant. The surfaces of the layer are at $z\!=\!\pm z_{{}_{\rm MAX}}$, and outside this ($|z|\!>\!z_{{}_{\rm MAX}}$) there is a hot gas with negligible density which exerts an external pressure $P_{{}_{\rm EXT}}$.

The layer is fully characterised by $P_{{}_{\rm EXT}}$, $c_{{}_{\rm S}}$ and its surface density $\Sigma_{{}_{\rm O}}$. A dimensionless parameter constructed from these quantities (Elmegreen & Elmegreen, 1978),

[TABLE]

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; $P_{{}_{\rm EXT}}\gg G\Sigma_{{}_{\rm O}}^{2}$) to 1 (self–gravity dominated).

From pressure equilibrium at the surfaces, $P_{{}_{\rm EXT}}\!=\!\rho(z_{{}_{\rm MAX}})c_{{}_{\rm S}}^{2}$, it follows

[TABLE]

Half–thickness of a stratified layer is defined by $H_{{}_{\rm HT}}\equiv\Sigma_{{}_{\rm O}}/(2\rho_{{}_{\rm O}})$. Equations (3), (1) and (2) then yield

[TABLE]

and the midplane density, $\rho_{{}_{\rm O}}$, depends on $\Sigma_{{}_{\rm O}}$ as

[TABLE]

Pressure dominated layers are of almost uniform density ($z_{{}_{\rm MAX}}\simeq H_{{}_{\rm HT}}$), while self–gravity dominated layers have a pronounced density maximum on the midplane ($z_{{}_{\rm MAX}}\gg H_{{}_{\rm HT}}$).

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 $\mn@boldsymbol{\omega}$ as a function of wavenumber $k$. In planar geometry, and for the dispersion relations in question, $\mn@boldsymbol{\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 $\mn@boldsymbol{\omega}$, which we denote for simplicity $\omega$, and we omit the oscillating imaginary part. The characteristic timescale of perturbation growth, its e–folding time, is defined as $t_{{}_{\rm EFOLD}}=1/\omega$.

For all the dispersion relations in question, there is a finite range of unstable wavenumbers $(0,k_{{}_{\rm MAX}})$, where $k_{{}_{\rm MAX}}$ is the highest unstable wavenumber. The maximum growth rate $\omega_{{}_{\rm FAST}}$ is attained for wavenumber $k_{{}_{\rm FAST}}$ which is always approximately $k_{{}_{\rm MAX}}/2$. The wavelengths corresponding to $k_{{}_{\rm MAX}}$ and $k_{{}_{\rm FAST}}$ are denoted $\lambda_{{}_{\rm MAX}}$ and $\lambda_{{}_{\rm FAST}}$, respectively. When describing a particular analytical estimate of the dispersion relation, the subscript begins with its name, e.g. $k_{{}_{\rm E78.FAST}}$ is the wavelength $k_{{}_{\rm FAST}}$ for E78. The abbreviations E78, V83 and W10 refer to the particular dispersion relation, not to the paper where they are firstly 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’s equation through the thickness of the shell. In planar geometry, the thin shell dispersion relation becomes

[TABLE]

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 $P_{{}_{\rm EXT}}$. Its semi-major axis is $r$, its semi-minor axis is the layer’s half-thickness, $H_{{}_{\rm HT}}$, and its total mass is $M$. Radial excursion of an element on the equator of the spheroid is regulated by the equation of motion

[TABLE]

(Boyd & Whitworth, 2005).

Wünsch et al. (2010) equate the instability growth rate to the radial contraction rate of the spheroid. Collapse from radius $r_{{}_{\rm O}}$, by a small factor $\epsilon$ to radius $(1-\epsilon)r_{{}_{\rm O}}$, in time $t_{\epsilon}$, fulfils $(\epsilon-1)r_{{}_{\rm O}}=\dot{r}_{{}_{\rm O}}t_{\epsilon}+\frac{1}{2}\ddot{r}_{{}_{\rm O}}t_{\epsilon}^{2}$. Using $\mn@boldsymbol{\omega}=1/t_{\epsilon}$ and $r=\lambda/2=\pi/k$, this yields

[TABLE]

The terms on the right hand side of Eq. (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 $\epsilon\!\sim\!0.1$, but, although $\epsilon$ affects the magnitude of $\mn@boldsymbol{\omega}_{{}_{\rm W10}}$, it does not influence the range of unstable wave-numbers.

2.2.3 Dispersion relation for a vertically stratified layer

Elmegreen & Elmegreen (1978) obtain the dispersion relation for a self-gravitating, pressure confined, semi-inifinite 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

[TABLE]

Here $c_{{}_{\rm EFF}}$ is the effective sound speed. $F_{{}_{\rm J}}$ and $F_{{}_{\rm D}}$ are reduction factors for self-gravity. Parameter $\eta$ 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 $\eta$ 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}}$, $\eta$, $F_{{}_{\rm J}}$ and $F_{{}_{\rm D}}$ are complicated functions of $k$, which can not generally be expressed in closed form.

The growth rate $\omega_{{}_{\rm E78.FAST}}$ of the most unstable wavenumber $k_{{}_{\rm E78.FAST}}$ in the limit $A\to 0$ (layers dominated by external pressure) is (eq. (47) in Kim et al. (2012))

[TABLE]

2.2.4 Comparison between the analytical estimates of the dispersion relation

The left panel of Figure 1 shows the dependence of $k_{{}_{\rm MAX}}$ normalised to $k_{{}_{\rm V83.MAX}}$ on parameter $A$ for different dispersion relations. For self–gravity dominated layers, both W10 and E78 predict almost identical ranges of unstable wavenumbers ($k_{{}_{\rm E78.MAX}}/k_{{}_{\rm W10.MAX}}\to 0.966$ as $A\to 1$), but V83 predicts a broader range with $k_{{}_{\rm W10.MAX}}/k_{{}_{\rm V83.MAX}}\to 0.518$. With increasing external pressure (decreasing $A$), the difference between E78 and V83 diminishes. However, when $A\la 0.5$, E78 extends to much higher wavenumbers, indicating that shorter wavelengths are unstable. And, as $A\to 0$, W10 and E78 have different limiting behaviour: $k_{{}_{\rm W10.MAX}}$ tends to a constant ($k_{{}_{\rm W10.MAX}}/k_{{}_{\rm V83.MAX}}\to 2.216$), whereas $k_{{}_{\rm E78.MAX}}\propto H_{{}_{\rm HT}}^{-1}\propto A^{-2}$.

Since the wavelength $\lambda_{{}_{\rm FAST}}$ has the highest growth rate, we assume that the fragment mass $M_{{}_{\rm FAST}}$ is approximately the mass confined inside a circle of radius $\lambda_{{}_{\rm FAST}}/2$, i.e. $M_{{}_{\rm FAST}}=\pi\Sigma_{{}_{\rm O}}(\lambda_{{}_{\rm FAST}}/2)^{2}=\pi\Sigma_{{}_{\rm O}}(\pi/k_{{}_{\rm FAST}})^{2}$. Masses $M_{{}_{\rm FAST}}$ normalised to the midplane Jeans mass $M_{{}_{\rm J}}$ are plotted in the right panel of Fig. 1. V83 and W10 predict that when formed, the fragments are already gravitationally unstable ($M_{{}_{\rm FAST}}/M_{{}_{\rm J}}\ga 1$) for any $A$. On the other hand, E78 predicts gravitationally unstable fragments only when $A\ga 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-parallelised 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) which 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. 2. 2.

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 $\Delta a^{\mathrm{max}}$ from Equation 9 of Salmon & Warren (1994), or they use the approximation by Springel (2005),

[TABLE]

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. 3. 3.

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, 2. 2.

the cell contains minimum of gravitational potential, 3. 3.

all the gas inside a sphere of accretion radius $r_{{}_{\rm ACC}}$ centered at the cell is above a specified density threshold, 4. 4.

the gas inside the sphere is gravitationally bound, Jeans-unstable and converging, 5. 5.

radius $r_{{}_{\rm ACC}}$ 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 ($H_{{}_{\rm HT}}$; see Eqn. 5 below), is always at least four grid cells (see Tables 1 to 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 (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\,\mathrm{M}_{{}_{\odot}}$ star, expanding into a medium with number density $n\!=\!10^{3}\,{\rm cm}^{-3}$. At time $t=0.81\,{\rm Myr}$, the shell has radius $R=3.86\,{\rm pc}$, and surface density $\Sigma_{{}_{\rm O}}=0.0068\,{\rm g}\,{\rm cm}^{-2}$. We use this value of $\Sigma_{{}_{\rm O}}$ for the initial conditions of all the simulations presented here, and vary $A$ by changing $P_{{}_{\rm EXT}}$. However, we note that the model is sufficiently simple that its physical parameters, $(\Sigma_{{}_{\rm O}},c_{{}_{\rm S}},P_{{}_{\rm EXT}})$, can be rescaled arbitrarily, as long as $A$ is unchanged, i.e. $P_{{}_{\rm EXT}}\propto\Sigma_{{}_{\rm O}}^{2}$.

The sound velocity, $c_{{}_{\rm S}}$, is related to the temperature $T$ by the ideal gas law $c_{{}_{\rm S}}^{2}\!=\!\gamma R_{gas}T/\mu$, where $\gamma$ is the effective barotropic exponent, $R_{gas}$ is the ideal gas constant, and $\mu$ is the mean molecular weight. Inside the layer, we set $\gamma=1.0001$ (effectively isothermal, $\gamma=1.0$ is excluded with the adopted numerical scheme) $\mu=2$ (pure molecular hydrogen), and $T=10\,{\rm K}$.

3.3 The external medium

We use two different kinds of the layer confinement: in Sections 4, 5 and 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 $\rho_{{}_{\rm AMB}}(z_{{}_{\rm MAX}})\lesssim 10^{-2}\rho_{{}_{\rm O}}$. In addition, to diminish the influence of sound waves reflection from borders of the computational domain, we extend the computational domain to $\sim\pm 3z_{{}_{\rm MAX}}$. Therefore, we set the ambient medium to be isothermal with temperature $T_{{}_{\rm AMB}}=300\,{\rm K}$ and mean molecular weight $\mu_{{}_{\rm 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.

[TABLE]

In the case of accreting layers, ram pressure is realised by supersonic accretion of gas with uniform density $\rho_{{}_{\rm ACC}}$, temperature $T_{{}_{\rm ACC}}$ and sound velocity $c_{{}_{\rm ACC}}$. The gas impacts at velocity $v_{{}_{\rm ACC}}$ (and hence Mach number $\mathscr{M}=v_{{}_{\rm ACC}}/c_{{}_{\rm S}}$) onto the upper surface of the layer at $z\simeq z_{{}_{\rm MAX}}$. The accreted gas is cold and molecular and the shock is isothermal, i.e. $T_{{}_{\rm ACC}}=T$ and $c_{{}_{\rm ACC}}=c_{{}_{\rm S}}$. The bottom surface at $z\simeq-z_{{}_{\rm MAX}}$ is confined by the ambient medium at temperature $T_{{}_{\rm AMB}}$, $T_{{}_{\rm AMB}}\gg 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. $\rho_{{}_{\rm ACC}}(v_{{}_{\rm ACC}}^{2}+c_{{}_{\rm ACC}}^{2})=\rho_{{}_{\rm 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 which 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 $\sigma$ calculated from the bottom side $z_{{}_{\rm BOT}}$ of the computational domain

[TABLE]

We cool to temperature $T$ any cell with the column density above threshold $\sigma_{{}_{\rm CRIT}}$. We set $\sigma_{{}_{\rm CRIT}}=2\int\limits_{z_{{}_{\rm BOT}}}^{-z_{{}_{\rm MAX}}}\,\rho_{{}_{\rm AMB}}(z^{\prime})\,dz^{\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 $\sigma(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 $T_{{}_{\rm AMB}}$ 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 $\Sigma$. 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 linearised equations, they are valid only as long as the perturbing amplitude $q_{1}$ of any quantity is smaller than its unperturbed value $q_{{}_{\rm O}}$. Accordingly, we define the linear regime of fragmentation if the maximum of the perturbed surface density $\Sigma_{1}$ is smaller than $\Sigma_{{}_{\rm 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{<}A{>}$$\_$${<}kH_{{}_{\rm HT}}{>}, 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 $kH_{{}_{\rm HT}}$. Thus, for example, simulation M18_020 treats a layer with $A=0.18$, and initial monochromatic perturbation $kH_{{}_{\rm HT}}=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 eqs. (20) – (23) in Kim et al. (2012)). The initial amplitude of perturbed surface density is $\Sigma_{1}(0)=0.01\Sigma_{{}_{\rm 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. Upper left panel of Figure 2 shows evolution of surface density perturbations for selected modes. Since a small perturbing amplitude in the surface density behaves as $\Sigma_{1}(t)=\Sigma_{1}(0)\mathrm{e}^{i\omega t}$, the instantaneous growth rate $\omega$ equals to the slope of the curves shown in the upper left panel of Fig. 2. $\omega$ is almost time independent throughout the simulations. Modes with $kH_{{}_{\rm HT}}\la 0.65$ are unstable and grow, modes with $kH_{{}_{\rm HT}}\ga 0.65$ are stable. We measure $\omega$ by linear fit to $\ln{(\Sigma_{1}/\Sigma_{1}(0))}$ over time interval $t^{fit}_{{}_{\rm O}}$, $t^{fit}_{1}$. The starting time $t^{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 timescale of sound crossing time through one half of the wavelength, i.e. $\simeq\pi H_{{}_{\rm HT}}/c_{{}_{\rm S}}$ (we use $t^{fit}_{{}_{\rm O}}=0.1\,\mathrm{Myr}$). The upper bound $t^{fit}_{1}$ is constrained by the condition for the perturbation to be small, i.e. $\Sigma_{1}<\Sigma_{{}_{\rm O}}$. Since model M18 terminates before this condition is fulfilled, we set $t^{fit}_{1}$ near the end of the simulation. We compare measured $\omega$ with the analytical dispersion relations Eq. (7), Eq. (9) and Eq. (10) in the upper right panel of Figure 2 and list measured e–folding time $t_{{}_{\rm NUM}}=1.0/\omega$ 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).

Monochromatic simulations for the self–gravity dominated layer ($A=0.99$, models M99) also show nearly time independent $\omega$ (lower left panel of Figure 2) with deviations towards the end only when the perturbing amplitude becomes large $\Sigma_{1}\ga\Sigma_{{}_{\rm O}}$. Measured $\omega$ ($t^{fit}_{{}_{\rm O}}=1.4\,\mathrm{Myr}$, $t^{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 timescales 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 $\|\vec{k}\|<k_{{}_{\rm O}}$ and are zero otherwise, so the amplitudes occupy a sphere in the Fourier space. To assure reasonable resolution, the highest wavenumber $k_{{}_{\rm O}}$ corresponds to the size of at least 4 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. $\mathrm{max}(k_{{}_{\rm E78.MAX}},k_{{}_{\rm V83.MAX}},k_{{}_{\rm W10.MAX}})<k_{{}_{\rm O}}$.

In order to compare fragmenting timescales between layers with different values of parameter $A$, it is necessary to choose comparable amplitudes of their initial perturbation. To reach the aim, we normalise all the perturbing amplitudes $\tilde{A}(k)$ to satisfy $\sqrt{\langle\tilde{A}(k)^{2}\rangle}/\tilde{A}(0)=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 normalisation constant. The lower limit is constrained by accuracy of the tree code. To fulfill these requirements, we set $\Sigma_{1}\simeq 0.1\Sigma_{{}_{\rm 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^{fit}_{{}_{\rm O}}$, $t^{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^{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^{fit}_{1}$ is constrained so as the modes with the highest growth rate have not reached the nonlinear regime yet. We find that $t^{fit}_{{}_{\rm O}}=1/\omega_{{}_{\rm E78.FAST}}$ and $t^{fit}_{1}=3/\omega_{{}_{\rm E78.FAST}}$ fulfill 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 panels of Figure 2) as for the monochromatic models presented in Section 4. The data are binned to reduce noise. The critical wavenumber $k_{{}_{\rm MAX}}$ as well as the growth rate of wavenumbers $k<k_{{}_{\rm MAX}}$ are again in a very good agreement with E78 both in self–gravity and external pressure dominated cases. Wavenumbers with $k>k_{{}_{\rm MAX}}$ are stable. We do not detect any significant deviation between our results and E78.

5.3 Evolution in non–linear regime

Evolution of surface density for self–gravity dominated model P99 is shown in Figure 3. The fragmentation begins with emergence of round objects (at time from 0.7Myr to 2.1Myr). 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.8Myr to 4.9Myr. Sink particles form in the densest parts of the filaments and accrete material from surrounding filaments. Comparing the last plot at 4.9Myr with a plot at the early stage of fragmentation (e.g. 2.1Myr), we see that when a clump appears, it monolithically collapses to a gravitationally unstable object.

Fragmentation in 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.3Myr to 1.0Myr) with masses smaller than the Jeans mass $M_{J}$. The subjeans masses are a direct consequence of fragmentation according to E78 for layers with low $A$ (right 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 equivalents 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.0Myr to 2.3Myr). Merging often leads to non–radial accretion resulting in spinning–up the fragments and disc formation around them (plot at time 2.3Myr). 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.

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 12. We experiment with several definitions of fragmenting time $t_{{}_{\rm FRG}}$, and conclude that taking $t_{{}_{\rm FRG}}$ 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 $t_{{}_{\rm FRG}}$ is also close to the time $t_{{}_{\rm SINK}}$ when sink particles in total contain 1/2 of the total mass of the layer (see Table 2 and left panel of Fig. 5).

Comparison between $t_{{}_{\rm FRG}}$ and E78 estimate $t_{{}_{\rm E78.FAST}}^{{}_{\rm EFOLD}}=1/\omega_{{}_{\rm E78.FAST}}$ (Eq. (11)) as a function of parameter $A$ is plotted in the left 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 self–gravity dominated case, but it underestimates the collapse time when external pressure dominates.

The parameters describing the layer are $\Sigma_{{}_{\rm O}}$, $c_{{}_{\rm S}}$ and $P_{{}_{\rm EXT}}$. From dimensional analysis it follows that any quantity with the dimension of time is given by

[TABLE]

where $f(A)$ is an unknown function of the parameter $A$. Since both $t_{{}_{\rm E78.FAST}}^{{}_{\rm EFOLD}}$ and $t_{{}_{\rm FRG}}$ are given by Eq. (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 $n_{{}_{\rm EFOLD}}(A)$ from our simulations, we can estimate $t_{{}_{\rm FRG}}$ for any layer because $t_{{}_{\rm E78.FAST}}^{{}_{\rm EFOLD}}$ is already known from Eq. (11).

The dependence of $n_{{}_{\rm EFOLD}}$ on $A$ is shown in the middle panel of Fig. 5 (solid line). Whereas $n_{{}_{\rm EFOLD}}$ 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 $t_{{}_{\rm NL}}$ (when $\Sigma_{1}=\Sigma_{{}_{\rm 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 $t_{{}_{\rm NL}}$ 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 $t_{{}_{\rm NL}}$, when $A$ is small, it takes many e–folding times to reach $t_{{}_{\rm FRG}}$. 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 panel of Figure 5. Analytical estimates for V83 and E78 ($M_{V83}$ and $M_{E78}$, respectively) and the Jeans mass $M_{{}_{\rm J}}$ for the midplane density (we use $M_{{}_{\rm J}}=4.26\;Ac_{s}^{4}/(G^{2}\Sigma_{{}_{\rm O}})$ from eq. (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=\pi\Sigma_{{}_{\rm O}}(\lambda_{{}_{\rm FAST}}/2)^{2}$ where $\lambda_{{}_{\rm FAST}}=2\pi/k_{{}_{\rm FAST}}$). 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\ga 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 subjeans 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 $\|\vec{k}\|\simeq k_{{}_{\rm FAST}}$ inclined at angles around $60\degr$ (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 idealised circumstances.

Name of the special models exploring possible pattern formation is in the form I{<}A{>}$$\_$${<}N_{w}{>}_[S|D], where the first two numbers after ”I” represent the parameter $A$ and the number after underscore the number $N_{w}$ 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 firstly 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 $\vec{k}_{i}$ ($i=1,2,3$) with $\|\vec{k}_{i}\|\simeq k_{{}_{\rm E78.FAST}}$ inclined at an angle around $60^{\circ}$ and amplify amplitudes of all modes inside a circle at centre $\vec{k}_{i}$ of radius $\mathrm{d}k$ by a bell shaped curve. Thus the perturbing amplitudes $\tilde{A}_{n}$ are calculated as

[TABLE]

where $\|\vec{k},\vec{k}_{i}\|_{2}=\sqrt{(k_{x}-k_{ix})^{2}+(k_{y}-k_{iy})^{2}}$ and $A_{int}$ and $\sigma_{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 $\mathrm{d}k$. We experimented with amplifying amplitudes of only the selected wavevectors $\vec{k}_{i}$ (model I99_3_D) and wavepackets inside non–zero radius $\mathrm{d}k$ centered on $\vec{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, $\vec{k}_{1}$ is amplified.

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\,\mathrm{Myr}$ and $2.90\,\mathrm{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(k_{x},k_{y})$. A regular pattern in the surface density would manifest itself as an anisotropy of function $B$. We look at the plane $(k_{x},k_{y})$ as it is given in polar coordinates specified by radius $\|\vec{k}\|$ and azimuthal angle $\alpha$. To search for the anisotropy, we analyse function $B$ evaluated on two different radially-averaged subsets in the $(k_{x},k_{y})$ plane. In the first case, $B$ is evaluated on a thin annulus with radius $\|\vec{k}\|=k_{{}_{\rm E78.FAST}}$, i.e. $B_{{}_{\rm RAD.1}}(\alpha)\equiv B(\|\vec{k}\|=k_{{}_{\rm E78.FAST}},\alpha)$. In the second case, the annulus is wider so it includes more modes around radius $k_{{}_{\rm E78.FAST}}$, i.e. $B_{{}_{\rm RAD.2}}(\alpha)\equiv B(k_{{}_{\rm E78.FAST}}/2\leq\|\vec{k}\|\leq 3k_{{}_{\rm E78.FAST}}/2,\alpha)$. If two modes inclined at an angle $\tilde{\alpha}$ have the highest amplitudes, function $B_{{}_{\rm RAD.J}}(\alpha)$ (letter $J$ stands for 1 or 2) has maxima separated by the angle $\tilde{\alpha}$. Thus the amplitudes of modes adjoined by the angle $\alpha$ are proportional to the value of azimuthal autocorrelation of $B_{{}_{\rm RAD.J}}(\alpha)$, i.e.

[TABLE]

Function $c_{J}(\alpha)$ detects interaction of modes inclined by any angle, not necessarily $60\degr$.

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 $c_{J}$. We denote $c_{{}_{\rm{MAX.}J}}$ the highest local maximum of function $c_{J}$ after its global maximum $c_{J}(0)$. During the simulations, the value of $c_{J}$ increases as amplitudes of modes increase. Since the interacting modes are predicted to grow faster than the rest of modes, we normalise $c_{{}_{\rm{MAX.}J}}$ with the mean $\langle c{{}_{J}}\rangle$.

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 Figure 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(k_{x},k_{y})$ and thus no formation of a regular pattern.

The bottom panel of Fig. 6 shows the angle $\alpha$ 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 on 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\degr$ 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 $a_{i}$ of the individual modes $\vec{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. eqs. (28) in (Fuchs, 1996)) with the initial mode amplitudes of model I99_3_D 222 Since the initial conditions for models I are not eigenfunctions, they do not exhibit the constant growth rate from the beginning. Eqs. (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 equal to that of the former at the end of the linear regime.

.The vertical bars mark the transition from the linear to non–linear regime. Recall that the growth rate $\omega$ 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.

In the previous paragraph, we demonstrate that the possible interaction between the triple of modes does not increase their growth rate above $\omega_{{}_{\rm 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, $\vec{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 $\alpha$ of the Fourier space. We do not investigate mode interaction in radial direction $\|\vec{k}\|$, 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{<}A{>}$$\_$${<}\mathscr{M}{>}$$\_{<}T_{{}_{\rm AMB}}{>}\_[N|$L]$. First number after letter A represents the value of the parameter $A$ at the beginning, followed by the Mach number $\mathscr{M}$ of the accreting medium and the temperature $T_{{}_{\rm AMB}}$ 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 $\mathscr{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 ($8\times$) smaller than in corresponding non–accreting model P30. Consequently, we can afford to use two times higher resolution $H_{{}_{\rm HT}}/\mathrm{d}z$, even though a number of grid cells in $x$ and $y$ directions are $4$ 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^{fit}_{{}_{\rm O}},t^{fit}_{1})=(1/\omega_{{}_{\rm E78.FAST}},2/\omega_{{}_{\rm E78.FAST}})$ for models with $T_{{}_{\rm AMB}}=300\mathrm{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 Figure 8. As the layer accretes, its surface density and parameter $A$ increase, we list these quantities at $t^{fit}_{{}_{\rm O}}$ and $t^{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^{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 $kH_{{}_{\rm HT}}$. The highest growth rate is around factor 2 higher than that of E78. To study possible influence of the temperature in the ambient medium, we calculate the same models with temperature $T_{{}_{\rm AMB}}=10^{4}\mathrm{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 $H_{{}_{\rm HT}}$; Table 4) than that for polychromatic models ($\simeq 4H_{{}_{\rm HT}}$; Table 2), which reproduce their analytic dispersion relation very well.

Evolution of density and velocity field inside the accreting layer (model A30_08_10_L) is shown in Figure 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 $kH_{{}_{\rm HT}}$ 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 toward 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 boundary conditions 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 ($\omega_{{}_{\rm E78.FAST}}t\simeq 4$), the growth of perturbations with $kH_{{}_{\rm HT}}\ga 0.7$ is saturated, so the range of unstable wavenumbers is comparable to that of the thermal–pressure confined 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 Figure 10. Model A30_05_03_N accretes at the highest rate (cf. the 6th column in Table 5) and substantial part of fragmentation occurs in self–gravity dominated case ($A\ga 0.6$; value of $A$ is at the upper right corner of each frame of Fig. 10). Fragments emerge at $t=1\rm{Myr}$ and subsequently collapse, forming filaments and then sink particles at their junctions ($t=1.5\rm{Myr}$ and $t=2\rm{Myr}$). Although the boundary conditions 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).

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 pressure–dominated case. At the beginning of the fragmentation ($t\la 1.0\rm{Myr}$), the density inside the layer is approximately constant and the overdensities seen in panels at $t=0.5\rm{Myr}$ and $t=1.0\rm{Myr}$ are corrugations of the contact discontinuity. The corrugations then gradually merge until they become gravitationally unstable and collapse (panel at $t=2.0\rm{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\,\mathrm{M}_{{}_{\odot}}$ to $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 pressure dominated case, its density (according to Eq. (6)) and therefore the Jeans mass $M_{{}_{\rm J}}$ is only a weak function of $A$. The highest change of $M_{{}_{\rm J}}$ is for the model with the highest accreting rate, i.e. A30_05_03_N; $M_{{}_{\rm J}}$ decreases only by factor $1.7\times$ 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 factor 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 boundary conditions. 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 which 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}$ ionising photons per second situated in a homogeneous medium with density $\rho_{{}_{\rm 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}\mathrm{s^{-1}}]$, $c_{s2}=[c_{{}_{\rm S}}/0.2\mathrm{km/s}]$, and $n_{3}=[n_{{}_{\rm ACC}}/10^{3}\mathrm{cm^{-3}}]$ where $c_{{}_{\rm S}}$ is sound velocity in the swept up shell, and $n_{{}_{\rm ACC}}$ is the number density in the surrounding medium. Following Whitworth et al. (1994b), we assume that fragmentation is accomplished in 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 factor 1.5 from their non–accreting counterpart (model P30). Mass of the fragments also differs at most by factor two 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 $\Sigma(t)$ and external pressure $P_{{}_{\rm EXT}}(t)$ fragments at any instant as the same layer with constant $\Sigma_{{}_{\rm O}}$ and $P_{{}_{\rm EXT}}$. Thus, when fragmentation starts at time $t_{{}_{\rm O}}$, it is accomplished at time $t_{1}$ fulfilling

[TABLE]

where the second equality follows from Eq. (11). Fitting measured data (Table 2 and middle panel in Fig. 5), we obtain

[TABLE]

We use our accreting models (with time dependent $\Sigma$) 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 Eq. (18) showing a good agreement.

The radius $R$ of a shell around an H ii region expands according to $R=R_{{}_{\rm ST}}(1+7c_{II}t/4R_{{}_{\rm ST}})^{4/7}$ (Spitzer, 1978), where $R_{{}_{\rm ST}}$ and $c_{II}$ is Strömgren radius and sound velocity in the ionised gas, respectively. Surface density of the shell and ram pressure due to swept up gas evolve as $\Sigma=\rho_{{}_{\rm ACC}}R(t)/3$ and $P_{{}_{\rm RAM}}=\rho_{{}_{\rm ACC}}(\mathrm{d}R/\mathrm{d}t)^{2}$. Thermal pressure $P_{{}_{\rm EXT}}$ inside the H ii region is equal to $P_{{}_{\rm RAM}}$. Thus approximating a small patch on the shell wall as a layer, its parameter $A$ (Eq. (3)) evolves as

[TABLE]

Evolution of $P_{{}_{\rm EXT}}$, $\Sigma$ and $A$ for walls of shells powered by sources with $\dot{N}_{49}=0.1$ (approximately an O9V star (Martins et al., 2005)), $\dot{N}_{49}=1$ (O6V star) and $\dot{N}_{49}=10$ (two O3V stars) expanding into a homogeneous ambient medium of particle density $n_{3}=0.1$, $n_{3}=1$ and $n_{3}=10$, respectively, is shown in Figure 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)\times(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 $t_{{}_{\rm O}}$ as determined from Eq. (22) below is plotted as the first asterisk on the evolutionary path. Following asterisks represent instants when integral in Eq. (18) reaches 0.25, 0.5, 0.75 and 1.0. They show that fragmentation proceeds in pressure dominated case.

Substituting $A$ from Eq. (20) for $t$ in Eq. (18) gives

[TABLE]

where $A_{{}_{\rm O}}=A(t_{{}_{\rm O}})$ and $A_{1}=A(t_{1})$ denote the value of parameter $A$ at the beginning and end of fragmentation, respectively.

We adopt $t_{{}_{\rm O}}$ 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\Sigma=R_{{}_{\rm J}}(4/(7t_{{}_{\rm O}}))^{2}$ (see appendix A in Whitworth et al. (1994b)) 333Our choice of $t_{{}_{\rm O}}$ is an approximation since for smaller fragments with $\lambda_{{}_{\rm MAX}}$, $\lambda_{{}_{\rm MAX}}<R_{{}_{\rm J}}$ self–gravity dominates stretching earlier than at $t_{{}_{\rm O}}$. However, the results (Eq. (23)) are weakly dependent on $t_{{}_{\rm O}}$; Any numerical coefficients in Eq. (23) would differ by less than 20 % if we used $t_{{}_{\rm O}}=0$ instead of Eq. (22). . We neglect the pressure gradient term since the unstable objects form via coalescence when the role of pressure support is already lost. Thus,

[TABLE]

8.2 The properties of fragments formed in a shell swept up by an expanding H ii region

Integrating Eq. (21) we find $A(t_{1})$ and from Eq. (20) fragmenting time $t_{1}$. We evaluate integral Eq. (21) for $c_{{}_{\rm S}}=0.2\mathrm{km/s}$ on a set of $\dot{N}_{49}\times n_{3}$ in range $(0.1,10)\times(0.1,10)$, and by polynomial fitting we find

[TABLE]

where $M_{{}_{\rm FRG}}$ is fragment mass, $r_{{}_{\rm FRG}}$ fragment radius, and $A_{{}_{\rm FRG}}$, $R_{{}_{\rm FRG}}$ and $\Sigma_{{}_{\rm FRG}}$ are parameter $A$, shell radius and surface density at time $t_{{}_{\rm FRG}}$, 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 Eq. (23b). We compare Eq. (23) with results of previous studies in Section 9.3.

Throughout this section, we have assumed that the shell fragments in external pressure dominated case. The highest value attained by $A$ in Eq. (23d) justifies this assumption ($A_{{}_{\rm FRG}}\la 0.6$) unless $n_{3}$ is either very high or $\dot{N}_{49}$ very low.

We restrict our study to sound velocity $c_{s2}=1$ (i.e. shell temperature $10\,\mathrm{K}$) to fulfill the condition $A_{{}_{\rm FRG}}\la 0.6$ for realistic values of $\dot{N}_{49}$ and $n_{3}$ because the area in $(\dot{N}_{49},n_{3})$ plane where the condition is met decreases as $c_{{}_{\rm S}}$ increases. Since fragmentation of self–gravity dominated shells is more complicated (time dependent $\rho_{0}$, gravitational focusing on emerging fragments, …), 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 $\Sigma_{1}$ is smaller than unperturbed $\Sigma_{{}_{\rm O}}$) clearly reproduce the semi–analytical estimate E78 in both 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 $k_{{}_{\rm MAX}}$ for $A\la 0.5$ (Fig. 1). V83 overestimates $k_{{}_{\rm MAX}}$ by a factor $2$ for $A=1$, predicts correct $k_{{}_{\rm MAX}}$ for $A\simeq 0.6$ and underestimates $k_{{}_{\rm MAX}}$ 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 $\Sigma_{{}_{\rm 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\la 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 $R_{{}_{C}}$ and central surface density $\Sigma_{1}$. The perturbed surface density decreases to zero from the centre to the edge at $R_{{}_{C}}$. 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 midplane. Then the condition for a lateral collapse is

[TABLE]

where $P_{1}$ and $P_{{}_{\rm O}}$ is the pressure in the centre and the edge of the perturbation, respectively. Factor $\beta$ takes into account that the perturbation generates non–spherically symmetric gravitational field. Factor $\beta$ is of order unity. Substituting $\Sigma_{1}+\Sigma_{{}_{\rm O}}$ and $\Sigma_{{}_{\rm O}}$ into Eq. (6) and using Eq. (5), Eq. (24) becomes $R_{{}_{C}}\beta\ga 2H_{{}_{\rm HT}}$. Assuming $R_{{}_{C}}=\lambda/2$, Eq. (24) takes the form

[TABLE]

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 (Eq. (9)) in $A\to 0$ predicts a significantly narrower range of unstable wavelengths, $kH_{{}_{\rm HT}}\leq 9\pi^{2}A^{2}/20$. The difference is due to how a perturbation is treated in the models. Imposing a perturbation with positive surface density $\Sigma_{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 midplane. The layer is almost of constant density where variations in surface density are dominated by corrugations of the surface. The weak dependence of $P_{1}$ on $\Sigma_{1}$ when $2P_{{}_{\rm EXT}}/(\pi G\Sigma_{{}_{\rm O}}^{2})\gg 1$ is apparent from Eq. (6). The density (pressure) gradient is much smaller than it would be if the perturbation were placed in the midplane, so the pressure gradient is easily overcame 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 $P_{{}_{\rm O}}$ and $P_{1}$, which needs a more massive fragment (with longer $\lambda_{{}_{\rm MAX}}$) to be overcame with its self–gravity.

In Section 8, assuming that shells swept up by H ii regions are at temperature $10{\rm K}$, we conclude that the shells fragment in 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 interstellar medium 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 $n_{{}_{\rm EFOLD}}$ on $A$ is done for the spectrum we adopt, the form of $n_{{}_{\rm EFOLD}}=n_{{}_{\rm EFOLD}}(A)$ is presumably a function of the spectrum shape. In other words, the spectrum introduces another dimensionless parameters into Eq. (15). The parameters may further alter the properties of the fragmenting shell (Eq. (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 timescale $t_{GC}$ proportional to its radius $r_{AC}$ and overdensity $\Sigma_{1}$ as $t_{GC}\simeq\sqrt{r_{AC}/(\pi G\Sigma_{1})}$. Given that the layer contains preexisting 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 center. 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 ionising 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 midplane; this assumption is violated in decelerating shells where it changes the dispersion relation as is shown by Iwasaki et al. (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 timescales. Taking into account approximations we made, we emphasize that Eq. (23) is rather an application of our results for layers than the definite answer to the problem of propagating star formation due to photoionising feedback. Assuming temperature of $10\mathrm{K}$ in the shells, Eq. (23) predicts that shell fragmentation tends to form objects of too low mass for star formation to propagate.

The dependence of our estimate Eq. (23) on $\dot{N}_{49}$ and $n_{3}$ is very close to the results provided by Whitworth et al. (1994a) and Iwasaki et al. (2011b) (with sound velocity $c_{s2}=1$). Compared to Whitworth et al. (1994a), Eq. (23a) predicts a slightly later fragmenting time ensuing in a larger shell radius (numerical constants they provide are $1.56\,\mathrm{Myr}$ and $5.8\,\mathrm{pc}$). However, our model predicts much smaller fragment masses ($3.1\,\mathrm{M}_{{}_{\odot}}$ vs. $23\,\mathrm{M}_{{}_{\odot}}$) and smaller fragment radii ($0.12\,\mathrm{pc}$ vs. $0.41\,\mathrm{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 et al. (2007) select one point in $\dot{N}_{49},n_{3}$ plane, $(\dot{N}_{49},n_{3})=(1.0,0.2)$ and follow its evolution by an SPH simulation. They find an approximate agreement with the work of Whitworth et al. (1994a). However, comparing density of their shell $\sim 10^{-21}\mathrm{g\;cm}^{-3}$ with the analytic value which is approx. $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 Eq. (23b) is close to the value $3.5\,\mathrm{M}_{{}_{\odot}}$ 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 case. 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\la 0.5$; Eq. (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; firstly it breaks into small gravitationally stable objects. It is an interesting feature because purely gravitational instability in planar geometry can form objects which are themselves stable against gravitational collapse. Then, the fragments continuously merge until they form an unstable object which 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 midplane of the layer.

In the non–linear regime, we investigate the possibility that the layer self–organises 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 factor $\simeq 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 ($10^{3}\,\mathrm{cm}^{-3}$) and temperature ($10\;\mathrm{K}$) in molecular clouds, the fragmentation is accomplished while the shell is still dominated by the external pressure. This leads to fragment masses $\simeq 3\,\mathrm{M}_{{}_{\odot}}$. 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 which 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“.

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 $\phi_{i}$ (Figure 12), we draw a sphere of the sink particle radius $r_{{}_{\rm SINK}}$ with centre at the potential minimum, and check whether the sphere contains at least one grid cell with negative binding energy $e_{{}_{\rm CELL}}$. The binding energy of a cell is

[TABLE]

where $m_{{}_{\rm CELL}}$, $\phi_{{}_{\rm CELL}}$, $v_{{}_{\rm CELL}}$ and $cs_{{}_{\rm CELL}}$ are mass, potential, bulk velocity and sound velocity of the cell. The kinetic energy is corrected by velocity of the mass centre $v_{{}_{\rm CL}}$ of the bound object, which is composed of all bound grid cells.

The gravitational binding energy is related to the potential at the saddle point $\phi^{ij}$ (the second letter in the superscript denotes the deeper minimum, i.e. $\phi_{j}<\phi_{i}$). The saddle for minimum $i$ is a pair of two neighbouring cells at the lowest potential $\phi_{{}_{\rm 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. 12). Mixed boundary conditions are taken into account, so the path can lead through a periodic boundary. No cell with $\phi_{{}_{\rm CELL}}>\phi^{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 Eq. (26) 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 firstly 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 can not be accessed later when investigating a deeper minimum.

We illustrate some of the features of the algorithm on Fig. 12. Note that the ordering of potential minima ($\phi_{m}>\phi_{l}>\phi_{i}>\phi_{j}>\phi_{k}$) is not reflected by ordering of saddles $\phi^{lk}>\phi^{jk}>\phi^{mj}>\phi^{ij}$. Only cells around $i$ with $\phi<\phi^{ij}$ can be assigned to minimum $i$. Since the saddle $\phi^{jk}$ for object $j$ is at a higher equipotential surface then previous saddle $\phi^{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$).