GRAVITATIONAL FRAGMENTATION OF EXPANDING SHELLS. II. THREE-DIMENSIONAL SIMULATIONS

We adopt the SPH method (e.g., Monaghan 1992). Since SPH is Lagrangian particle method, it is one of the best methods for problems where a wide low-density region and geometrically thin shell coexist. Instead of additional viscosity to handle shocks, we use the "Godunov SPH" where the results of the Riemann problem are used to calculate the interaction between SPH particles (Inutsuka 2002). The tree method (Barnes & Hut 1986) is used to calculate self-gravitational force.

2.1.1. Treatment of Contact Discontinuity

In the standard SPH method (Monaghan 1992), the ith particle has mass m_i, and its density is given by

$$\begin{equation} \rho _i = \sum _k m_k W({\bf x}_i-{\bf x}_k,\bar{h}_{ik}), \end{equation} \tag{ 1 }$$

where W(x, h) is the kernel function and we use the Gaussian kernel, h is the smoothing length, and $\bar{h}_{ik}$ is an average between h_i and h_k. The equation of motion of the ith particle is given by

$$\begin{equation} \frac{\Delta {\bf v}_i}{\Delta t} = - \sum _k m_k \left(\frac{P_i}{\rho _i^2} + \frac{P_k}{\rho _k^2} \right) \frac{\partial }{\partial {\bf x}_i}W({\bf x}_i - {\bf x}_k, \bar{h}_{ik}). \end{equation} \tag{ 2 }$$

The standard SPH has a shortcoming when calculating the CD between a hot and a cold gas (Agertz et al. 2007). Agertz et al. (2007) reported that the artificial tension stabilizes the Kelvin–Helmholtz instability with different densities in the standard SPH. In contrast, Cha et al. (2010) have shown that Godunov SPH can describe the Kelvin–Helmholtz instability reasonably well.

In addition, Tartakovsky & Meakin (2005) and Hu & Adams (2006) proposed a new method that treats the CD naturally in the context of a multi-phase fluid. They modified Equations (1) and (2) to treat the CD as follows:

$$\begin{equation} \rho _i = m_i \sum _k W({\bf x}_i-{\bf x}_k,\bar{h}_{ik}) \end{equation} \tag{ 3 }$$

and

$$\begin{eqnarray} \frac{\Delta {\bf v}_i}{\Delta t} &=& {-} \frac{1}{m_i}\sum _k \left[ P_i\left(\frac{m_i}{\rho _i}\right)^2 + P_k\left(\frac{m_k}{\rho _k} \right)^2\right]\nonumber \\ &&\times \frac{\partial }{\partial {\bf x}_i} W({\bf x}_i - {\bf x}_k. \bar{h}_{ik}). \end{eqnarray} \tag{ 4 }$$

The term m_i/ρ_i represents the volume of the ith particle ∼h³_i. Therefore, if we determine the particle mass of hot and cold gases so that its smoothing length distributes smoothly across the CD, the terms inside the bracket in Equation (4) is smooth across the CD, leading to much better behavior at the CD. Their new method is very useful in problems where the position of the CD is known in advance as in our modeling. Equation (4) satisfies the relation of m_iΔv_i = −m_kΔv_k, suggesting that the momentum conservation is guaranteed. Furthermore, we apply the Godunov technique to Equation (4) as

$$\begin{eqnarray} \frac{\Delta {\bf v}_i}{\Delta t} &=& {-} \frac{1}{m_i}\sum _k P^*_{ik}\left[ \left(\frac{m_i}{\rho _i}\right)^2 + \left(\frac{m_k}{\rho _k} \right)^2\right]\nonumber \\ &&\times \frac{\partial }{\partial {\bf x}_i} W({\bf x}_i - {\bf x}_k, \bar{h}_{ik}), \end{eqnarray} \tag{ 5 }$$

where P*_ik is the result of the Riemann problem between the ith and the kth particles (Inutsuka 2002).

Massive stars emit ultraviolet photons (hν>13.6 eV) and produce H ii regions around them. Here, we consider a massive star that emits ionizing photons with photon number luminosity Q_UV(s⁻¹) into the ambient gas with uniform density of ρ_E = mn_E, where n_E and m are the number density and the mean mass of the ambient gas particle, respectively.

We construct as simple a model as possible to concentrate on the physics of GI of expanding shells. Figure 1 illustrates the schematic picture of our model. In this paper, we do not calculate the radiative transfer of ionizing photons but introduce hot and cold gases that are assumed to evolve keeping their temperatures constant. The hot gas is occupied in r < R_CD (the dotted region), and the thin shell is in R_CD < r < R_SF, where R_CD and R_SF are positions of the CD and the SF, respectively. The pre-shock ambient cold gas in r>R_SF is assumed to be uniform. In our model, the inner boundary is set at r = R_b inside the hot bubble (R_b < R_CD). In the H ii region, the detailed balance between the recombination and the ionization is approximately established all the time. Therefore, the interior pressure of the H ii region is given by

$$\begin{equation} P_\mathrm{int} = \rho _\mathrm{E} c_\mathrm{II}^2 \left(\frac{R_\mathrm{ST}}{R_\mathrm{CD}} \right)^{3/2}, \end{equation} \tag{ 6 }$$

where R_ST is the Strömgren radius,

$$\begin{equation} R_\mathrm{ST} = \left(\frac{3Q_\mathrm{UV}}{4\pi \alpha _\mathrm{B} n_\mathrm{E}^2} \right)^{1/3}. \end{equation} \tag{ 7 }$$

We assume that the pressure of the H ii region is spatially constant because of the high temperature. The hot gas at r = R_b is pushed by the interior pressure P_int.

Zoom In Zoom Out Reset image size

If we simulate the overall shell, the required number of SPH particles, N_tot, is too large to calculate. Therefore, to save N_tot, we calculate a part of the shell. The schematic picture of the calculation domain is shown in Figure 2 and is designated by |tan⁻¹(y/x)| ⩽ θ_b and |tan⁻¹(z/x)| ⩽ θ_b. Since the solid angle subtended by the calculation domain from the center O is given by

$$\begin{equation} \Omega _\mathrm{b} \simeq 4 \theta _\mathrm{b}^2\;\;\;\mathrm{for}\;\;\theta _\mathrm{b}\ll 1, \end{equation} \tag{ 8 }$$

the total number of SPH particles can be reduced by a factor of Ω_b/4π ∼ θ²_b/π compared with the calculation of the whole shell.

Zoom In Zoom Out Reset image size

In this paper, for convenience, the units of the time, length, and mass scales are taken to be

$$\begin{equation} t_0 = \sqrt{\frac{3\pi }{32G\rho _\mathrm{E}}} = 1.6\;n_\mathrm{E,3}^{-1/2}\;\mathrm{Myr}, \end{equation} \tag{ 9 }$$

$$\begin{equation} R_0 = \left(\frac{7 c_\mathrm{II} t_0}{\sqrt{12}} \right)^{4/7} R_\mathrm{ST}^{3/7} = 5.9\; Q_\mathrm{UV,49}^{1/7}\: n_\mathrm{E,3}^{-4/7}\;\mathrm{pc}, \end{equation} \tag{ 10 }$$

and

$$\begin{equation} M_0 = \rho _\mathrm{E} R_0^3 = 5.0\times 10^3\; Q_\mathrm{UV,49}^{3/7}\: n_\mathrm{E,3}^{-5/7}\,M_\odot, \end{equation} \tag{ 11 }$$

respectively, where Q_UV,49 = Q_UV/10⁴⁹ s⁻¹ and n_E,3 = n_E/10³ cm^-3. The dependence of the expansion law on the parameters (Q_UV,  n_E) are approximately eliminated by using the non-dimensional quantities normalized by t₀, R₀, and M₀ (see Figure 1 in Paper I).

It is difficult to resolve a very thin shell in early evolutionary phase even if we only calculate part of the shell (see Figure 2). Thus, we use a grid-based 1D hydrodynamical code to describe the early evolutionary phase of the H ii region expansion, and we switch to 3D SPH at t = 0.4 t₀. As a grid-based numerical method, we use the 1D spherically symmetric Lagrangian Godunov method (van Leer 1997). The innermost mesh at R_b is pushed by the interior pressure given by Equation (6). When the 1D simulation reaches t/t₀ = 0.4, the radial profiles of physical quantities are used as the initial condition of the 3D calculations. We smooth the distribution of the sound speed across the CD in the scale of the smoothing length. It is assumed that the temperature of the hot gas is 64 times higher than that of the cold gas. The specific value of the hot gas temperature does not change our conclusion as long as it is much larger than the temperature of the cold gas.

We calculate the expansions of H ii regions around the 41 M_☉ (the high-mass (HM) model) and 19 M_☉ (the low-mass (LM) model) stars that are embedded by the uniform ambient gas of n_E = 10³ cm⁻³. Simulation parameters are listed in Table 1. The temperature of the cold gas T_c is assumed to be 10 K. The opening angles of simulation domains are set to θ_b = 2π/26 and 2π/14 for the HM and the LM models, respectively, so that the calculation domains contain a single wavelength of the most unstable mode predicted from the thin-shell linear analysis by Elmegreen (1994). The mass of SPH particles are set so that the initial thickness is resolved by five SPH particles. Since the thickness increases with time, the relative resolution becomes better as the shell expands. In the later phase, the thickness is resolved by ∼15 particles. The total number of SPH particles for the HM and the LM models are 4.00 × 10⁶ and 2.26 × 10⁶, respectively.

Corrugation-type perturbations are put into the shell at the initial state. We displace the SPH particles so that their densities do not change, or ∇ · Δx_i = 0, where Δx_i is the displacement of the ith particle. The functional form of the displacement of the shell is assumed to be ∝−cos(lϕ), where ϕ = tan⁻¹(y/x). Therefore, the shell has a negative displacement at ϕ = 0. The SPH particles are displaced, keeping their velocity, the sound speed, and the mass fixed. We concentrate on the evolution of a single mode in each calculation. Dependence on l is investigated by many runs.

The moving boundary condition at r = R_b(t) is realized by using "ghost particles" (Takeda et al. 1994). The time evolution of R_b(t) is obtained with the results of the 1D simulations. At the four surface areas (y = ±xtan θ_b, z = ±xtan θ_b) in the pyramid-shaped calculation domain (see Figure 2), the rotational periodic boundary condition is imposed.

Since the gravitational force is a long-range force, we have to take into account the particles in the whole solid angle. However, we only have information in the part of the solid angle as shown in Section 2.3. In this paper, the gravitational force is evaluated using the following method. Figure 3 shows the schematic picture of the gas sphere viewed from the x-direction for the case with θ_b = 2π/10. The latitude and longitude lines are plotted at intervals of θ_b. The region is surrounded by the thick solid lines representing the calculation domain (see Figure 2). We rotate the calculation domain iθ_b degrees in the θ-direction, where i = −π/(2θ_b), ..., π/(2θ_b). At each rotation step in the θ-direction, we rotate the domain kθ_b degrees in the ϕ-direction, where k = 1, ..., 2π/θ_b. After that, the information in the whole solid angle can be obtained. However, for i ≠ 0, the rotated domains in the ϕ-direction have overlaps. This comes from the fact that the pyramid cannot fill 3D space. For example, the dashed regions in Figure 3 correspond to the rotated domains in the θ-direction. We remove the overlapped SPH particles with |tan⁻¹(y/x)|>θ_b from the rotated domains. In the simulation, we rotate not the SPH particles but the tree structure, and calculate the gravitational force from the rotated tree structures at each rotation step.

Zoom In Zoom Out Reset image size

With this method, the supposed particle distribution in the whole solid angular 4π is not completely cyclic. However, the effect of the approximation is negligible since θ_b is very small in this paper. For i = ±1 (near the equatorial plane), the number of removed SPH particles is negligible compared with the total number. On the other hand, for large |i| (near the poles), the fraction of removed SPH particles becomes large. However, since the rotated domain is very far from the calculation domain, the detailed particle distribution is not important.

In this paper, we neglect the gravitational force outside the SF to prevent the ambient gas from collapsing towards the center. We find the radius R_grav of the most distant particle from the center that has density ρ>ρ_th, where ρ_th is a threshold density (Bisbas et al. 2009). In our simulations, we adopt ρ_th = 1.1ρ_E. Only SPH particles within R_grav are assumed to feel the gravitational force that is calculated with the above method. This means that there is a discontinuity in the gravitational force at R_grav. In the Appendix, we investigate the effect of the discontinuity and confirm that it does not influence our results as long as R_grav is inside the shock transition layer. The main reason is that the pressure suddenly changes in the shock transition layer, suggesting that the pressure gradient is much larger than the jump in the gravitational force.

In order to evaluate the growth of perturbations in the SPH calculations, we introduce indicators for fluctuations of density and the position of the CD. The solid angle Ω_b is divided into N_Ω cells, and the direction of the ith cell is defined by the unit vector n_i (1 ⩽ i ⩽ N_Ω). The radial profiles of the physical quantities along n_i are obtained from the SPH calculations. Here, we define the angle-dependent maximum density ρ_max(n_i) along n_i. The position of the CD along n_i is determined as the position where the temperature coincides with a critical value that is set to $\sqrt{2}T_\mathrm{c}$. The specific choice of the critical value is not important in our results. We average ρ_max(n_i) and R_CD(n_i) over all directions. The mean value of Q = (ρ_max, R_CD) is defined by

$$\begin{equation} \langle Q \rangle = \frac{1}{N_\Omega }\sum _{i} Q({\bf n}_i). \end{equation} \tag{ 12 }$$

As the indicators of the fluctuations, we evaluate the dispersions of these quantities:

$$\begin{equation} \Delta Q \equiv \sqrt{\frac{1}{N_\Omega }\sum _{i} \left\lbrace Q({\bf n}_i) - \langle Q\rangle \right\rbrace ^2}. \end{equation} \tag{ 13 }$$

As a test calculation, we present the evolution of the shell without initial perturbations in the HM model. Figure 4 shows the time evolution of 〈R_CD〉 (Equation (12)) evaluated from the SPH simulation (the open circles). For comparison, we plot the result of the 1D simulation with the solid line. One can see that the SPH simulation can reproduce the results of the 1D calculation very well.

Zoom In Zoom Out Reset image size

Figure 5 shows snapshots of density (the upper panel) and pressure (the lower panel) profiles for t/t₀ = 0.5, 0.7, 1.0, and 1.26. The abscissae indicate the distance from the CD at each time. Only SPH particles in |y|/R₀ < 0.05 and |z|/R₀ < 0.05 are plotted by the dots. For comparison, the density profiles in the 1D simulation are superimposed with the thick gray lines in the upper panel. It is seen that the result of the SPH calculation agrees with that of the 1D simulation very well. Our 3D simulation clearly represents the structure and evolution of the shell although it is very thin.

Zoom In Zoom Out Reset image size

Note that the pressure profiles in the lower panel of Figure 5 are smooth at the CD (r − 〈R_CD〉 = 0) thanks to the modified Godunov SPH shown in Section 2.1.1. In the standard SPH, a wiggle in the pressure profile appears at the CD. In order to see whether the pressure profile is smooth even when perturbation is added, the density and pressure profiles for t = 1.0 and l = 52 are plotted in Figure 6. The detailed investigation of the results with perturbation is shown in Section 4. The abscissa is the distance from 〈R_CD〉. We plot SPH particles in ϕ₀ − Δϕ < ϕ < ϕ₀ + Δϕ and |z|/R₀ < 0.05, where Δϕ = 0.001θ_b. Figures 6(a) and (b) correspond to ϕ = 0 where δρ has the maximum value and ϕ₀ = θ_b/2 where δρ has the minimum value, respectively. One can clearly see that the pressure profile is smooth at the CD even with perturbation.

Zoom In Zoom Out Reset image size

The density profile of the decelerating shell is asymmetric with respect to the density peak. Iwasaki et al. (2011) developed a semi-analytic method for deriving the time evolution of the density profile by assuming that the shell reaches hydrostatic equilibrium along the pressure gradient, the inertia force owing to the deceleration, and the radial gravitational force toward the center (Whitworth & Francis 2002). Iwasaki et al. (2011) found that the time evolution of the density profile can be divided into three evolutionary phases: deceleration-dominated, intermediate, and self-gravity-dominated phases (see Figure 2 in Paper I). In the early phase (deceleration-dominated phase), the density peak is at the SF because the inertia force is larger than the gravitational force. As the shell expands, the inertia force decreases while the gravitational force increases. Thus, the density peak moves from the SF toward the CD as shown in Figure 5. In the intermediate phase, the density peak is closer to the SF than the CD. In the later phase (self-gravity-dominated phase), the density peak is closer to the CD than the SF.

It is also found that the density profiles for various sets of (n_E,  Q_UV) are characterized by a single parameter, that is, the typical Mach number,

$$\begin{equation} {\cal M}_0 = \frac{4}{7}\frac{R_0}{c_\mathrm{s} t_0}= 7\; Q_\mathrm{UV,49}^{1/7}\: T_\mathrm{c,10}^{-1/2}\: n_\mathrm{E,3}^{-1/14}, \end{equation} \tag{ 14 }$$

where T_c,10 = T_c/10 K.

They found the development of the GI is strongly influenced by two effects, asymmetry of the density profile and boundary conditions.

The asymmetry of the density profile strongly influences the GI, especially, in the self-gravity-dominated phase (for example, the shell at t/t₀ = 1.26 in Figure 5). In this later phase, the distance from the density peak to the SF is larger than H₀, where $H_0=c_\mathrm{s}/\sqrt{2\pi G\rho _{00}}$ is the scale height of the shell and ρ₀₀ is the peak density. On the other hand, the distance from the density peak to the CD is smaller than H₀. The gas tends to collect toward the density peak because the unperturbed gravitational potential has the minimum value there. The gas near the SF can collapse toward the density peak because the sound wave cannot travel from the density peak to the SF. This mode is the "compressible mode." On the other hand, the gas near the CD cannot collapse to the peak. However, the GI can proceed even there through the deformation of the CD that makes the gravitational potential deeper. This mode is the "incompressible mode" (Elmegreen & Elmegreen 1978; Lubow & Pringle 1993). Therefore, the GI of the shell in the later phase has properties of compressible and incompressible modes.

Moreover, Iwasaki et al. (2011) found that the growth rate is enhanced compared to the symmetric case through the combination of the GI and the Rayleigh–Taylor instability.

The expanding shell is confined by the SF on the leading surface and by the CD on the trailing surface. Dispersion relations of Elmegreen (1994, E94) and Whitworth et al. (1994) are essentially based on the dispersion relation of the layer confined by the SF on both surfaces. The shock-confined layer is stabilized by the tangential flow generated by the obliqueness of the SF compared with the layer confined by the thermal pressure of the hot gas. Therefore, the linear analysis in Paper I that takes into account SF+CD boundary conditions gives larger growth rate compared with E94.

In summary, as described in Sections 3.2 and 3.3, Paper I predicts that the GI begins to grow earlier than predicted in E94. The development of the GI in the later phase is accompanied by the significant deformation of the CD.

First, we present the results of the HM model with perturbations. We calculate four runs for the case with l = 26, 52, 78, and 156. Figures 7(a), (b), (c), and (d) represent the time evolution of the perturbation amplitude for the wavenumbers l = 26, 52, 78, and 156, respectively. The thick solid and the thick dashed lines correspond to the density perturbation Δρ_max/〈ρ_max〉 and the deformation of the CD 30 × ΔR_CD, respectively. In each panel, Δρ_max/〈ρ_max〉 and ΔR_CD grow in a similar way. Perturbations with l = 52–78 grow with the largest growth rate. For larger angular wavenumber l = 156 and smaller angular wavenumber l = 26, they grow more slowly. The prediction from E94 is also superimposed by the thin dashed line in each panel. It is found that in our 3D results, perturbations begin to grow earlier than predicted by E94.

Zoom In Zoom Out Reset image size

Let us compare the results of SPH simulations with the linear analysis in Paper I. However, Paper I neglected the evolutionary effects, i.e., expansion of the shell and the mass accretion through the SF. To include these effects approximately, we modify the instantaneous growth rate from Γ to Γ_evo as follows:

$$\begin{equation} \Gamma _\mathrm{evo}(l,t) \equiv -3\frac{V_\mathrm{s}(t)}{R_\mathrm{s}(t)} + \sqrt{\left(\frac{V_\mathrm{s}(t)}{R_\mathrm{s}(t)} \right)^2 + \Gamma (l,t)^2}, \end{equation} \tag{ 15 }$$

where V_s is the velocity of the shell. This modification is inspired by E94. The time evolution of R_s(t) and V_s(t) is determined by the thin-shell model shown in Section 2 in Paper I. One can see that the V_s/R_s terms correspond to the evolutionary effect of the shell that stabilizes the GI. The prediction based on Paper I is superimposed by the thin solid line in each panel of Figure 7. It is evaluated by exp [∫_tΓ_evo(l, t')dt']. From Figure 7, one can see that the linear analysis of Paper I describes the growth of perturbations very well.

Figure 8 shows the time sequence of cross sections through the z = 0 plane for l = 52. The color scale indicates density, and the arrows represent velocity relative to the fluid with the maximum density. The abscissa and the ordinate axes correspond to the azimuthal angle ϕ, and the distance from the CD divided by 〈R_CD〉, respectively. Figure 8(a) represents the initial condition where the shell around ϕ = 0 is displaced to the negative radial direction. Figure 8(b) shows that the tangential flow is generated by the obliqueness of the SF, and collects the gas around ϕ = 0. As a result, enhanced pressure pushes the SF outward in the radial direction (see Figure 8(c)). The gravitational potential around ϕ = 0 becomes deep, leading to gas accumulation (see Figure 8(d)). This mode of the GI is similar to the incompressible mode in the sense that the deformation develops at the boundary (Elmegreen & Elmegreen 1978; Lubow & Pringle 1993). In the later phase (Figure 8(e)), the CD highly deforms while the SF is almost flat. This is consistent with prediction of Paper I (see Section 3). Figure 8(e) is very similar to Figure 11 in Paper I that shows the eigenfunction in the parameters of the HM model at t/t₀ = 1.3 for l = 52.

Zoom In Zoom Out Reset image size

Figure 9 is the same as Figure 8 but for a larger wavenumber case, l = 156. The initial state is shown in Figure 9(a). Because of large wavenumber, the tangential flow collects the gas more quickly in Figure 9(b), which is similar to Figure 8(c). In this moment, since the SF at ϕ = 0 deforms to the positive radial direction, the gas moves away from ϕ = 0 due to the tangential flow. As a result, at t/t₀ = 0.8 (Figure 9(c)), the SF and the CD are almost flat while the velocity field exists inside the shell. The density distribution at t/t₀ = 1.0 (Figure 9(d)) has the opposite phase compared to the initial state. In this moment, the perturbation begins to grow but the growth rate is smaller than that for l = 52. The CD deforms largely at the final state in both Figures 8 and 9. In this sense, the behavior of the GI is similar to that for l = 52.

Zoom In Zoom Out Reset image size

Next, we present results of the LM model. Because the central star is less energetic, the Mach number of the expanding shell is smaller than that in the HM model, suggesting that the shell is geometrically thicker. Therefore, perturbations are expected to grow more slowly than the HM model. This feature is confirmed in Figure 10, which is similar to Figure 7. The most unstable mode is found in a smaller wavenumber, l ∼ 28, and the growth rate is smaller than in the HM model. We found that perturbations begin to grow earlier than predicted by E94 also in this case. The linear analysis in Paper I describes the GI very well. Figure 11 shows the cross section through z = 0 for l = 56 at t/t₀ = 1.5. The significant deformation of the CD is seen as in the HM model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We presented only two examples of the HM and the LM models in Section 4. In this section, we predict the GI in the shells with various parameters (n_E, Q_UV, T_c) by using the results of the linear analysis presented in Paper I. As shown in Figures 7 and 10, it is found that the results of the SPH simulations are well described by the growth rate, Γ_evo. Using this growth rate, one can predict the time t_bgn when the GI begins to grow using the condition Γ_evo(l_bgn,  t_bgn) = 0, where l_bgn is the angular wavenumber of the most unstable mode at t = t_bgn. This condition is rewritten as

$$\begin{equation} \sqrt{8}\frac{V_\mathrm{s}}{R_\mathrm{s}} = \Gamma (l_\mathrm{bgn},\:t_\mathrm{bgn}), \end{equation} \tag{ 16 }$$

where the left-hand side corresponds to the evolutionary rate owing to expansion and mass accretion, and the right-hand side corresponds to the growth rate of the GI. At the early phase, the evolutionary rate is larger than the growth rate of the GI, indicating that the shell is stable. As the shell expands, V_s/R_s decreases while Γ increases. Therefore, the GI begins to grow when the growth rate of the GI is larger than the evolutionary rate. The radius of the shell at t = t_bgn is defined by R_s,bgn.

Figures 12(a), (b), (c), and (d) show contour maps of t_bgn, R_s,bgn, λ_max ≡ 2πR_s,bgn/l_bgn, and M_bgn ≡ Σ(t = t_bgn)π(λ_bgn/2)² in the (n_E,  Q_UV) plane derived by using linear analysis, respectively. Here, T_c is assumed to be 10 K. From Figure 12, one can see that t_bgn, R_s,bgn, λ_bgn, and M_bgn decrease with increasing n_E. This can be understood by scaling the typical number density of the shell, $n_\mathrm{sh}=n_\mathrm{E} {\cal M}_0^2$, that corresponds to the shocked gas density. From Equation (14), we have n_sh ∝ n^6/7_E. Therefore, for larger n_E, the shell becomes denser and unstable earlier (Figure 12(a)), and the most unstable scale is smaller (Figures 12(c) and (d)). Next, let us see the dependence of t_bgn, R_s,bgn, λ_bgn, and M_bgn on Q_UV. Since the central star is more energetic for larger Q_UV, the Mach number of the expanding shell $({\cal M}_0\propto Q_\mathrm{UV}^{1/7})$ is larger, indicating that the shell is denser and is more unstable. Therefore, t_bgn, λ_bgn, and M_bgn decreases with increasing Q_UV. On the other hand, R_s,bgn increases with Q_UV because the expansion velocity increases with Q_UV.

Zoom In Zoom Out Reset image size

From Figure 12, since the contour lines are almost straight in the plane (log₁₀n_E,  log₁₀Q_UV), the quantities are expected to have power-law dependence on Q_UV and n_E, or equivalently, on ${\cal M}_0$. In Figure 13, we plot t_bgn/t₀, R_s,bgn/R₀, and l_bgn as a function of ${\cal M}_0$ for various parameters that correspond to the open circles in Figure 12. The open circles correspond to the results of the linear analysis. From Figure 13, one can see that the results of the linear analysis are well described by the analytic formulae, $t_\mathrm{bgn}/t_0=1.62{\cal M}_0^{-0.75}$, $R_\mathrm{s,bgn}/R_0=1.2{\cal M}_0^{-0.375}$, and $l_\mathrm{bgn}=3.2{\cal M}_0^{1.5}$ which are plotted by the solid lines. Using Equation (14), the analytic formulae are rewritten as

$$\begin{equation} t_\mathrm{bgn} = 0.6\:Q_\mathrm{UV,49}^{-0.11}\: T_\mathrm{c,10}^{0.375}\: n_\mathrm{E,3}^{-0.45}\;\mathrm{Myr}, \end{equation} \tag{ 17 }$$

$$\begin{equation} R_\mathrm{s,bgn} = 3.4\: Q_\mathrm{UV,49}^{0.09}\: T_\mathrm{c,10}^{0.19}\: n_\mathrm{E,3}^{-0.54}\;\mathrm{pc}, \end{equation} \tag{ 18 }$$

and

$$\begin{equation} l_\mathrm{bgn} = 54\; Q_\mathrm{UV,49}^{0.21}\: T_\mathrm{c,10}^{-0.75}\: n_\mathrm{E,3}^{-0.11}. \end{equation} \tag{ 19 }$$

Zoom In Zoom Out Reset image size

The wavelength λ_bgn and the mass M_bgn that correspond to l_bgn are given by

$$\begin{equation} \lambda _\mathrm{bgn} = 0.4\; Q_\mathrm{UV,49}^{-0.125}\: T_\mathrm{c,10}^{0.94}\: n_\mathrm{E,3}^{-0.44}\;\mathrm{pc}, \end{equation} \tag{ 20 }$$

and

$$\begin{equation} M_\mathrm{bgn} = 3.5\; Q_\mathrm{UV,49}^{-0.16}\: T_\mathrm{c,10}^{2.06}\: n_\mathrm{E,3}^{-0.42}\,{M_\odot }. \end{equation} \tag{ 21 }$$

Whitworth et al. (1994, W94) also derived similar formulae by using the thin-shell linear analysis similar to E94. The dependence of Equations (17)–(21) on parameters Q_UV, T_c, and n_E shows close agreement to that in W94 although the numerical factors are different. Our results show that the GI begins to grow earlier and as a result, the fragment mass is smaller than their results. In detail, t_bgn, R_bgn, and λ_bgn are roughly half of those in W94, and M_bgn is roughly 1/8 of that in W94. This is because the stabilization effect of the SF does not work on the trailing surface that is the CD as shown in Paper I. As pointed out in W94, the properties of fragments are insensitive to Q_UV, and mainly depend on n_E and T_c. Indeed, in Equations (20) and (21), the dependence of λ_bgn and M_bgn on n_E are close to the Jeans length and the Jeans mass of the shell (λ_J,  M_J) ∝ n^−1/2_sh ∝ n^−3/7_E.

Elmegreen (1989) investigated the GI in a decelerating shocked layer. He numerically integrated perturbation equations with respect to time by taking into account the time evolution of the layer and by averaging physical quantities over the thickness. He has taken into account boundary conditions on the SF and CD. He called the boundary effect on the CD "pinching force." In his linear analysis, the destabilizing effect of the pinching force does not work efficiently. One reason is that he used Equation (3) in Paper I as the pressure of the H ii region though the geometry is plane-parallel. Therefore, the pressure at the CD is underestimated, leading us to underestimate the pinching force of the CD. The other reason is the geometry of the gravitational potential. The asymmetry of the density profile of the spherical shell is qualitatively different from that of the plane-parallel layer. In the plane-parallel geometry, the static symmetric layer peaks mid-plane. If the layer decelerates, the density peak is always closer to the SF than the CD. Therefore, the tangential flow behind the SF erodes the density perturbation more efficiently. On the other hand, for the case with shells, the peaks can be closer to the CD than the SF (see Section 3.1). From these reasons, he derived lower growth rates than us.

Dale et al. (2009) simulated the gravitational fragmentation of expanding shells confined by temporary constant thermal pressure of hot rarefied gases on both sides. In their calculation, the motion of the shell is determined so that the pressure at the leading surface is the same as that at the trailing surface. Therefore, the density peak is always around the mid-plane of the shell, and the density profile is almost symmetric. In their calculation, the column density decreases with time because the shell expands, keeping the mass fixed. Therefore, the pressures at the boundaries approach the peak pressure. They found that the confining pressure accelerates fragmentation in the later phase, and describe this effect as "pressure-assisted" gravitational fragmentation. This mode is the same as the incompressible mode. Recently, taking into account the effect of the thickness of the shell approximately, Wünsch et al. (2010) established a semi-analytic linear analysis that explains results of their SPH simulations. Different from them, in Paper I, we have taken into account the effect of SF+CD boundary condition approximately.

We discuss other potentially important effects that are not analyzed in this paper.

5.3.1. Vishniac Instability

5.3.2. Thermal Instability

In this paper, we did not solve the radiative transfer equation of the ionizing photons. Instead, we assumed that the trailing surface is the CD instead of the ionization front (IF). There are two effects of the ionization on the GI. One is that a fraction of the shell becomes ionized gas as the H ii region expands. Therefore, our simulations overestimate the column density of the shell. The swept-up mass by the SF is given by

$$\begin{equation} M_\mathrm{sw} = \frac{4\pi }{3}\rho _\mathrm{E} R_\mathrm{s}^3. \end{equation} \tag{ 22 }$$

On the other hand, the ionized mass is given by

$$\begin{equation} M_\mathrm{ion} = \frac{4\pi }{3}\rho _\mathrm{i} R_\mathrm{s}^3, \end{equation} \tag{ 23 }$$

where ρ_i is the density of the ionized gas, and we neglect the difference between the radii of the IF and SF. The real shell mass is given by M_sw − M_ion. Inside the H ii region, the balance between ionization and recombination is approximately established. Therefore, the density of the H ii region becomes

$$\begin{equation} \rho _i = \left(\frac{3m^2 Q_\mathrm{UV}}{4\pi \alpha _\mathrm{B} R_\mathrm{s}^3} \right)^{1/2} = \rho _\mathrm{E}\left(\frac{R_\mathrm{ST}}{R_\mathrm{s}} \right)^{3/2}, \end{equation} \tag{ 24 }$$

where m is the mean weight of the gas particles and we use Equation (7). Therefore, from Equations (22)–(24), the ratio of M_ion to M_sw is

$$\begin{equation} {\cal R}\equiv \frac{M_\mathrm{ion}}{M_\mathrm{sw}} = \left(\frac{R_\mathrm{ST}}{R_\mathrm{s}} \right)^{3/2}. \end{equation} \tag{ 25 }$$

If ${\cal R}\ll 1$, the effect of the ionization is negligible. In Figure 1 in Paper I (see also Figure 4), the Strömgren radius is as large as 0.1R₀ in various parameters (Q_UV, n_E). At the starting time of our simulations, t = 0.4 t₀ (R_s ∼ 0.6R₀), ${\cal R}\sim (0.1/0.6)^{3/2}=0.07$. As the shell expands, ${\cal R}$ decreases. At t = t₀, ${\cal R}\sim (0.1/10)^{3/2}\sim 0.03$. Therefore, the contribution of the ionization is limited to only several percent.

The other is that the IF induces an instability on the shell. The D-type IF is analogous to a combustion front that is unstable to corrugational deformations (Landau & Lifshitz 1987). Vandervoort (1962) found that the D-type IF is unstable. In this instability, the most unstable scale is comparable to the thickness of the IF. Moreover, Garcia-Segura & Franco (1995) showed that the VI is strongly modified by the IF. The shell rapidly fragments and finger-like structures are generated in their two-dimensional simulations. Note, however, that they assumed the power-law density distribution ∝r^−w of the ambient gas with the uniform density core in the center with 0.2 pc, where the parameter w is set to 2. The expansion of the IF depends on the power-law index w in the 1D model. If w>1.5, even if the SF emerges in front of the IF in the uniform density core, the IF quickly gets ahead of the SF, and eventually the whole of the shell and the ambient gas are ionized. Moreover, Hosokawa (2007) investigated the expansion of the IF and the dissociation front, and showed that star formation is expected to be suppressed when w>1. This is because the column density decreases as the shell expands and the FUV photon easily escapes in front of the shock. Therefore, the model by Garcia-Segura & Franco (1995) corresponds to the case where the triggered star formation is not simply expected. We need to investigate the expansion of the IF for the case with the shallower density profile w < 1 in detail using the radiative multi-dimensional simulations taking into account self-gravity. Moreover, Williams (1999) discovered shadowing instability in the R-type IF before the emergence of the SF. This instability may also disturb the shell in the very early phase of evolution.

It is well known that the magnetic field is important in the dynamics of the ISM although it is neglected in this paper. Nagai et al. (1998) investigated the GI of a pressure-confined isothermal layer threaded by uniform magnetic field by using linear analysis. They found that the magnetic field cannot stabilize the GI. Moreover, they found the layer fragment in filamentary clouds whose longitudinal axis is parallel or perpendicular to the magnetic field depending on the thickness of the layer. However, their analysis is restricted to the static magnetized layer. It is important to investigate the stability of the magnetized shocked layer or shell in dynamical modeling.