A&A 374, 746-755 (2001)
DOI: 10.1051/0004-6361:20010768

Gravitational instability of expanding shells

Solution with nonlinear terms

R. Wünsch - J. Palouš

Astronomical Institute, Academy of Sciences of the Czech Republic, Bocní II 1401, 141 31 Praha 4, Czech Republic

Received 9 October 2000 / Accepted 22 May 2001

Abstract
A model of the thin shell expanding into a uniform ambient medium is developed. Density perturbations are described using equations with linear and quadratic terms, and the linear and the nonlinear solutions are compared. We follow the time evolution of the fragmentation process and separate the well defined fragments. Their mass spectrum is compared to observations and we also estimate their formation time.

Key words: hydrodynamics - instabilities - stars: formation - ISM: bubbles - ISM: kinematics and dynamics - galaxies: ISM

1 Introduction

H$_{\alpha}$ emission, which is frequently observed along the periphery of giant expanding HI shells, is related to star formation, which has been triggered due to the compression of the ISM. As an example of correlated star formation, the ring of young stars observed near to the Sun consisting of Orion, Perseus, Sco-Cen and other OB associations can be used (Lindblad et al. 1997). The HI hole No. 35 in the dwarf galaxy IC 2574 (Walter & Brinks 1999) with the H$_{\alpha}$ emission found along its rims serves as another example of triggered star formation.

The analysis of the stability of a pressure-confined slab performed by Elmegreen & Elmegreen (1978) has been extended to spherical expanding shocks by Vishniac (1983). The time development of the gravitational collapse of linear perturbations in decelerating, isothermal shocked layers has been examined numerically and analytically by Elmegreen (1989, 1994). The purely hydrodynamical nonlinear instability has been discussed by Vishniac (1994). In this paper, we continue with nonlinear analysis of the gravitational instability of spherically symmetric shells expanding into stationary homogeneous medium.

We modify the approach adopted by Fuchs (1996) who described the fragmentation of uniformly rotating self-gravitating disks. If some conditions are fulfilled, the expanding shell may become gravitationally unstable and break to fragments. The inclusion of higher order terms helps to determine with better accuracy than the linear analysis when, where and how quickly it happens.

The Rayleigh-Taylor (R-T) instability is not expected to develop in the situation explored because a spherical shock expanding into the homogeneous interstellar medium is always decelerated collecting stationary ambient medium. R-T instabilities may be important in different situations when the density of the ambient medium drops down sufficiently quickly so that the shell can accelerate mixing hot and cold gas components. This is the case of very active SF regions where the shell interacts with previously formed fragments. The region 30 Doradus in LMC may serve as an example of R-T instability in action as described by Redman et al. (1999).

For highly supersonic flows multiple shocks may develop (Falle 1981), but at that time the shells are gravitationally stable due to squeezing connected to the fast expansion. Later, when they decelerate to velocities less than 50 kms^-1, the gravity starts to be important, while radiative instabilities of the outer shock described by Strickland & Blondin (1995) loose their influence.

This work may be extended to nonspherical oscillations using the formalism worked out by Bicák & Schmidt (1999) for cosmological applications. Here we also ignore the deviations from spherical symmetry resulting from initial asymmetry of the energy input. In a smooth medium with only large scale density gradients the shell approaches quickly the spherical symmetry as demonstrated by Bisnovatyi-Kogan & Blinnikov (1982). Nonradial perturbations resulting from inhomogeneities of the ambient medium and variations of the shell surface density will be discussed in a subsequent paper.

2 The expanding shell in a static, homogeneous medium

The energy input from an OB association or other sources creates a blast-wave propagating into the ambient medium (Ostriker & McKee 1988; Bisnovatyi-Kogan & Silich 1995). Since its radius is during the majority of the evolution much larger than its thickness, the thin shell approximation can be used. The blast-wave is considered as an expanding infinitesimally thin layer surrounding the hot medium inside. The analytical solution of the expansion in the thin shell approximation was derived by Sedov (1959). In a static, homogeneous medium without the external or internal gravitational field, the blast-wave is always spherically symmetric, it sweeps the ambient matter and decelerates. Its evolution is given by the equation of motion

$${{\rm d} \over {\rm d}t} ( M V ) = S ( P_{\rm int} - P_{\rm ext} ),$$ (1)

the equation of mass conservation

$${{\rm d} \over {\rm d}t} M = S \rho_0 V,$$ (2)

and the equation of state

$$P_{\rm int} = {2 \over 3} {E_{\rm th} \over V {\it ol}},$$ (3)

where M and V are the total mass collected in the shell and its expansion velocity, S is the shell surface, $P_{\rm int}$ and $P_{\rm ext}$are pressures inside and outside of the shell, $\rho_0$ is the density of the ambient medium and $E_{\rm th}$ and $V {\it ol}$ are the thermal energy and volume inside of the shell. Neglecting the external pressure and keeping the $E_{\rm th} = {\rm const.}$, and assuming $E_{\rm th} = {3 \over 5} E_{\rm tot}$, where $E_{\rm tot}$ is the total input energy from the source, the radius of the shell R grows with time as

 

$$\left({R \over {\rm pc}}\right)\ = 72.2 \left(\frac{ E_{\rm tot}}... ...o_0}{ 1 {\rm cm}^{-3}}\right)^{-1/5} \left(\frac{ t}{ {\rm Myr}}\right) ^{2/5},$$ (4)

it decelerates as

 

$$\left(\frac{ V}{ {\rm km\,s}^{-1}}\right)\ = \ 28.2 \left(\frac{ ... ...o_0}{ 1 {\rm cm}^{-3}}\right)^{-1/5} \left(\frac{ t}{ {\rm Myr}}\right)^{-3/5},$$ (5)

and its unperturbed surface density $\Sigma_0$ grows with time as

 

$$\left({\Sigma_0 \over {\rm {\it M}_{\odot}}\, {\rm pc}^{-2}}\righ... ...\rho_0}{ 1 {\rm cm}^{-3}}\right)^{4/5} \left(\frac{t}{ {\rm Myr}}\right)^{2/5}.$$ (6)

The linear analysis of the gravitational instability of expanding shells by Elmegreen (1994) gives for the instantaneous maximum growth rate $\omega_{\rm BGE}$of a transversal perturbation of a shell

$$\omega_{\rm BGE} = -{3 V \over R} + \left[{V^2 \over R^2} + \left( \pi G \Sigma_0 \over c \right)^2\right]^{1/2},$$ (7)

where c is the sound speed within the shell. The instability occurs for $\omega_{\rm BGE} > 0$. Inserting Eqs. (4-6) into Eq. (7), we may derive the time $t_{\rm b}$ when the instability occurs for the first time. $\omega _{\rm BGE}(t_{\rm b}) = 0$ for

$$\left({t_{\rm b} \over {\rm Myr}}\right) = 28.8 \left(\frac{ c}{ ... ... {\rm erg}}\right)^{-1/7} \left(\frac{ \rho_0}{ 1 {\rm cm}^{-3}}\right)^{-4/7}.$$ (8)

The ratio of the wavelength $\lambda$ of the fastest perturbation to R

$$\lambda / R = {2 c^2 \over {G R} \Sigma_0}$$ (9)

is at $t \ge t_{\rm b}$ less than ${\pi c \over \sqrt{2} V}$. The sound speed cin the dense and cold shell is always smaller than the speed of sound in the ambient medium, which is smaller than the expansion speed of the shell V. Therefore, $\lambda / R << 1$.

In a subsequent paper, we shall also discuss the instability in non-spherical shells: the values of $R,\ V,$ and $\Sigma_0$ will be taken from numerical simulations.

3 Hydrodynamical and Poisson equations on the surface of the shell

We consider the cold and thin shell of radius R surrounding the hot interior and expanding with velocity V into a uniform medium of density $\rho_0$. The intrinsic surface density of the shell $\Sigma$ is composed of unperturbed part $\Sigma_0$ plus the perturbation $\Sigma_1$( $\Sigma=\Sigma_0 + \Sigma_1$). Perturbation $\Sigma_1$ results from the flows on the surface of the shell redistributing the accumulated mass. We assume that $\Sigma_0$ corresponds to R as $\Sigma_0= \rho_0 R/3$, which means that all the encountered mass is accumulated to the shell. (It comes from $\Sigma_0 = \frac{4} {3}\pi R^3 \rho_0 / 4\pi R^2$.)

The mass conservation law in a small area on the surface of the shell is

$$\frac{\partial m}{\partial t}+\left( \nabla, m\vec v \right) = A\ V \rho_0,$$ (10)

where m is mass in the area $A = 4\pi\alpha R^2$, $\alpha$ is a fixed small fraction of the sphere. $\vec v$ denotes a two-dimensional velocity of surface flows connected to perturbations above the stretching of the area due to expansion of the unperturbed shell: $\vec r$ and $\vec v = \dot{\vec r}$ are two-dimensional vectors in the tangential plane of the shell at the central point of the considered area. We consider angular coordinates $\vec\Theta=\vec r/R$ and angular velocity $\vec\Omega = \vec v/R$ to describe the evolution on the surface of the shell (see Fig. 1).

	Figure 1: The coordinates on the shell surface: $\vec\Theta=\vec r/R$for the position and the angular velocity $\vec\Omega = \vec v/R$for the surface flows.
Open with DEXTER

With $\Sigma = m / A$ we obtain continuity equation in a form

$$\frac{\partial \Sigma}{\partial t}+2\Sigma\frac{V}{R} +\Sigm... ... \right)+R\left(\vec\Omega, \nabla\right) \Sigma = V \rho_0 .$$ (11)

The equation of motion for flows on the shell surface has a form

$$\frac{1}{A}\frac{{\rm d}(m\vec v)}{{\rm d}t}=-c^2\nabla\Sigma-\Sigma\nabla \Phi,$$ (12)

where c is the constant isothermal sound speed inside the cold shell, $\Phi$is the gravitational potential generated by the mass distribution given in the tangential plane to the shell. The shell is confined by the pressure of the bubble interior and the outer shock. We assume that the pressure distribution through the shell corresponds to constant isothermal sound speed c and constant volume density inside of the shell with a slight deviations close to its inner edge, which is in contact with the hot interior. Using the continuity Eq. (11), we can rewrite the Eq. (12) as

 

$$R\frac{\partial\vec\Omega}{\partial t} + V\vec\Omega+3V\vec\Omega... ...a, \nabla \right) \vec\Omega - R^2 \vec\Omega \left( \nabla, \vec\Omega \right)$$

$$=-\frac{c^2}{\Sigma}\nabla\Sigma-\nabla\Phi.$$ (13)

The gravitational potential $\Phi$ is related to the surface density by the Poisson equation

$$\Delta\Phi = 4\pi G\Sigma\delta(z)\ ,$$ (14)

where G is the constant of gravitation and $\delta(z)$ is a delta function of the space coordinate z perpendicular to the surface of the shell.

4 Perturbation analysis

We assume a small perturbation of the shell surface density $\Sigma_1 \ll \Sigma_0$ which evolves due to surface flows given with velocity $\vec v$. The perturbed hydrodynamical Eqs. (11), (13) and the perturbed Poisson Eq. (14) have form

$${\partial \Sigma_1 \over \partial t} + 2\Sigma_1 {V \over R}... ... \right) + R \left( \vec \Omega, \nabla \right) \Sigma_1 = 0,$$ (15)

 

$$R {\partial \vec \Omega \over \partial t} + R^2 \left( \vec\Omega, \nabla \right) \vec\Omega - R^2 \vec\Omega \left( \nabla, \vec\Omega \right)$$

$$= -{c^2 \over \Sigma_0}(1-{\Sigma_1 \over \Sigma_0}) \nabla \Sigma_1 - \nabla \Phi_1 - 4V\vec \Omega + 3V \vec \Omega {\Sigma_1 \over \Sigma_0},$$ (16)

$$\Delta \Phi_1 = 4\pi G \Sigma_1 \delta (z),$$ (17)

where the $1/\Sigma$ in Eq. (13) was evaluated up to quadratic terms in $\Sigma_1$: ${c^2 \over \Sigma}\nabla \Sigma = {c^2 \over \Sigma_0} (1 - {\Sigma_1 \over \Sigma_0} +$ higher order terms) $\nabla \Sigma_1$, and $\nabla \Sigma_0 = 0$ on the shell surface.

The perturbation of the surface density $\Sigma_1$ and the angular velocity of the surface flows $\vec\Omega$ can be written as

 

$$\Sigma_1 \Sigma_{10} + \sum_{\vec\eta} \Sigma_{\vec\eta} {\rm e}^{i \left(\vec{\eta}, \vec{\Theta} \right) } \ ,$$

$$\vec\Omega \vec\Omega_0+\sum_{\vec\eta} \vec\Omega_{\vec\eta} {\rm e}^{i \left( \vec{\eta}, \vec{\Theta}\right)} \ ,$$ (18)

which may be inserted to the perturbed Eqs. (15) and (16). $\vec\eta$ denotes a dimensionless wave-vector $\vec\eta=\vec kR$. We assume no surface macroscopic flow through all the considered area which means $\vec\Omega_0=0$. Further we assume $\Sigma_{10}=0$ (mass accumulation due to expansion to the ambient medium is included in $\Sigma_0$). The Fourier transform of the Eq. (15) is

 

$$\dot\Sigma_{\vec\eta} + \Sigma_0\left(-i \vec\eta, \vec\Omega_{-\... ...}\left(\vec\Omega_{\vec\eta-\vec\eta'}, -i\vec\eta'\right) \Sigma_{\vec\eta'}=0$$ (19)

where we used the identity $\vec\eta'+(\vec\eta-\vec\eta')=\vec\eta$.

The solution of the Poisson Eq. (17)

$$\nabla\Phi_{1} = -2\pi G\sum_{\vec\eta}\Sigma_{\vec\eta} \fra... ...ert\vec\eta\vert} {\rm e}^{i\left(\vec\eta, \vec\Theta\right)}$$ (20)

can be inserted to the Fourier transform of the Eq. (16). We get:

 

$$R\dot{\vec\Omega}_{\vec\eta} + R\sum_{\vec\eta'} \left( \vec\Omeg... ...ght) \vec\Omega_{\vec\eta'} - R\sum_{\vec\eta'} \vec\Omega_{\vec\eta-\vec\eta'}$$

$$\left(-i\vec\eta', \vec\Omega_{-\vec\eta'} \right) = - \frac{c^2}... ...ma_0^2}\sum_{\vec\eta'}i\vec\eta' \Sigma_{\vec\eta-\vec\eta'}\Sigma_{\vec\eta'}$$

$$+ 2\pi G\frac{i\vec\eta}{\vert\vec\eta\vert}\Sigma_{\vec\eta} - 4... ...V}{\Sigma_0}\sum_{\vec\eta'}\vec\Omega_{\vec\eta'} \Sigma_{\vec\eta-\vec\eta'}.$$ (21)

5 Linear analysis

We give the solution of Eqs. (19) and (21) using the linear terms only. This linear solution will be later compared to the results obtained with the nonlinear terms. The complete linear analysis of the expanding shell was also done by Elmegreen (1994). Linearized Eqs. (19) and (21) have a form

$$\dot\Sigma_{\vec\eta} + \Sigma_0 \left( -i\vec\eta, \vec\Omega_{-\vec\eta}\right) + 2\frac{V}{R} \Sigma_{\vec\eta} = 0$$ (22)

$$R\dot{\vec\Omega}_{\vec\eta} = - \frac{c^2}{\Sigma_0}i\vec k ... ...ert\vec\eta\vert}\Sigma_{\vec\eta} - 4V\vec\Omega_{\vec\eta}.$$ (23)

Angular velocity $\vec\Omega_{\vec\eta}$ can be split in two components parallel and orthogonal to the wave-vector $\vec\eta$

$$\vec\Omega_{\vec\eta}=\Omega_{\vec\eta_\parallel} \frac{\vec... ...\eta_\perp}\frac{\vec\eta_\perp}{\vert\vec\eta_\perp\vert} \ ,$$ (24)

where $\vec\eta_\perp$ is a vector in the tangential plane, which is perpendicular to the wave-vector $\vec\eta$and $\vert\vec\eta_{\perp}\vert=\vert\vec\eta\vert$.

We get the set of equations

$$\dot{\Sigma}_{\vec\eta}=-2\frac{V}{R}\Sigma_{\vec\eta} -i\eta\Sigma_0\Omega_{\vec\eta_\parallel}$$ (25)

$$\dot{\Omega}_{\vec\eta_\parallel}=(i\frac{2\pi G}{R} -i\frac... ...^2})\Sigma_{\vec\eta} -4\frac{V}{R}\Omega_{\vec\eta_\parallel}$$ (26)

$$\dot{\Omega}_{\vec\eta_\perp}=-4\frac{V}{R}\Omega_{\vec\eta_\perp}$$ (27)

which can be formally written as

$$\pmatrix{ \dot\Sigma_{\vec\eta} \cr \dot\Omega_{\vec\eta_\p... ...cr \Omega_{\vec\eta_\parallel} \cr \Omega_{\vec\eta_\perp} }.$$ (28)

We find the eigenvalues and eigenvectors of the linear operator $\cal L$.

$$\omega_{\eta}^{(1,2)} = i\frac{3V}{R} \pm \sqrt{-\frac{V^2}{R^2} + \frac{\eta^2 c^2}{R^2} - \frac{2\pi G\Sigma_0\eta}{R}}$$ (29)

$$\omega_{\eta}^{(3)}=i4\frac{V}{R} \ ,$$ (30)

where $\eta$ is the magnitude of $\vec\eta$.

The related eigenvectors are

$$\pmatrix{\Sigma_{\vec\eta} \cr \Omega_{\vec\eta_\parallel} \... ...ta\Sigma_0 \cr i 2\frac{V}{R} - \omega^{(1,2)}_{\eta } \cr 0}$$ (31)

$$\pmatrix{\Sigma_{\vec\eta} \cr \Omega_{\vec\eta_\parallel} \cr \Omega_{\vec\eta_\perp}}^{(3)} = \pmatrix{0\cr 0\cr 1} .$$ (32)

The eigenvalues $\omega_{\eta}$ are time dependent and can be used to obtain a criteria for the instability of the shell. For a short time for which the change of the linear operator $\cal L$ is small (and we can always find the time interval which is short enough), the solution of Eqs. (28) is a part of the exponential function $\sim$ ${\rm e}^{i\omega t}$. The $\omega^{(3)}_{\eta}$ has always meaning of the decrease of perturbations. If $\omega^{(1,2)}_{\eta}$ have got a real part, solution is stable with decreasing oscillations. If not, $\omega^{(1)}_{\eta}$ indicates decrease, $\omega^{(2)}_{\eta}$ can be imaginary negative and it have meaning of the perturbations growth rate. The time evolution of the imaginary part of the $\omega^{(2)}_{\eta}$ is shown by Fig. 2.

	Figure 2: The time dependence of the imaginary part of the $\omega ^{(2)}$, which can cause the instability. The Sedov solution was used with following parameters: total energy $E_{\rm tot} = 10^{53}$ erg, density of ambient medium $n_0 = 1~{\rm cm}^{-3}$, average molecular weight $\mu = 1.3$, sound speed in the shell $c = 1~{\rm km\,s}^{-1}$.
Open with DEXTER

Since the eigenvalues $\omega_{\eta}^{(1,2)}$ depend on $\eta$ through the relation (29), the maximum perturbation growth rate $\omega_{\eta, {\max}}^{(1,2)}$ can be found:

$$\omega_{\eta,\max}^{(1,2)} = i\frac{3V}{R} \pm \sqrt{-\frac{V^2}{R^2} - \frac{\pi^2 G^2 \Sigma_0^2 }{c^2}}$$ (33)

and occurs at the wavenumber:

$$\eta_{\max}=\frac{\pi G \Sigma_0 R}{c^2}\cdot$$ (34)

$\omega^{(2)}_{\eta, \max}$ is close to $\omega_{\rm BGE}$ derived by Elmegreen (1994), actually $\omega^{(2)}_{\eta, \max} = - i \omega_{\rm BGE}$. If the shell is unstable, i.e. the imaginary part of the $\omega^{(2)}_{\eta}$ is negative, for the certain time, the shell may break to fragments. In Ehlerová et al. (1997) is the fragmentation integral for the maximum perturbation growth rate defined as:

$$I_{\rm f}(t) \equiv \int_{t_{\rm b}}^{t} i\omega^{(2)}_{\eta,{\rm max}}(t') {\rm d}t' \ ,$$ (35)

where $t_{\rm b}$ is the time when the instability begins. We generalize it for any unstable mode with the wavenumber $\eta$ using an analogy to the value for $\omega_{\rm BGE}$ as

$$I_{{\rm f}, \eta }(t) \equiv \int_{t_{\rm b}}^{t} i\omega^{(2)}_{\eta }(t') {\rm d}t'.$$ (36)

The fragmentation time $t_{{\rm f}, \eta }$(the time giving some development level of a fragment) is defined as the time when the fragmentation integral is equal to one:

$$I_{{\rm f}, \eta }(t_{{\rm f}, \eta }) = \int_{t_{\rm b}}^{t_{{\rm f}, \eta }} i\omega^{(2)}_{\eta } (t') {\rm d}t' = 1\ .$$ (37)

6 Nonlinear analysis

The equations with nonlinear terms are solved by the similar procedure as adopted by Fuchs (1996). We rewrite the nonlinear Eqs. (19) and (21) in the form

$$\pmatrix{ \dot\Sigma_{\vec\eta} \cr \dot\Omega_{\vec\eta_\p... ...ec\eta_\parallel} \cr \Omega_{\vec\eta_\perp} } + {\cal N} \ ,$$ (38)

where $\cal L$ is the linear part and ${\cal N}$ represents the non-linear terms. We search for a solution of Eqs. (38) as a combination of the eigenvectors obtained from the previous linear analysis

 

$$\pmatrix{ \Sigma_{\vec\eta} \cr \Omega_{\vec\eta_\parallel} \cr \... ..._{\vec\eta} \cr \Omega_{\vec\eta_\parallel} \cr \Omega_{\vec\eta_\perp}}^{(3)},$$ (39)

where $\psi_{\vec\eta(t)}$, $\xi_{\vec\eta(t)}$ and $\phi_{\vec\eta(t)}$ are time dependent amplitudes of eigenvectors. There are always four solutions differing in eigenvectors multiplied by $\pm i$ or $\pm 1$. We select only solutions with physical relevance given by eigenvectors (31). We find orthonormal vectors $(\Sigma_{\vec\eta},\Omega_{\vec\eta_\parallel}, \Omega_{\vec\eta_\perp})^{(1,2,3)}$in order that

$$(\Sigma_{\vec\eta}, \Omega_{\vec\eta_\parallel},\Omega_{\vec\eta_... ...\vec\eta} \cr \Omega_{\vec\eta_\parallel} \cr \Omega_{\vec\eta_\perp} }^{(j)}=1$$

$$(\Sigma_{\vec\eta}, \Omega_{\vec\eta_\parallel},\Omega_{\vec\eta_... ...a} \cr \Omega_{\vec\eta_\parallel} \cr \Omega_{\vec\eta_\perp} }^{(k\not=j)}=0,$$ (40)

where j, k = 1,2,3. The orthonormal vectors are

 

$$\begin{array}{l} (\Sigma_{\vec\eta},\Omega_{\vec\eta_\parallel},\... ...ta_\parallel},\Omega_{\vec\eta_\perp})^{(3)}= \left( 0,0,1 \right) \end{array}.$$ (41)

We insert ansatz (39) into Eq. (38), multiply it by the orthonormal vectors (41) and obtain a set of equations for amplitudes $\psi_{\vec\eta(t)}$, $\xi_{\vec\eta(t)}$ and $\phi_{\vec\eta(t)}$

 

$$\dot\psi_{\vec\eta} i\omega^{(1)}\psi_{\vec\eta} - \left( \frac{\partial}{\partial t}... ...\vec\eta_\parallel},\Omega_{\vec\eta_\perp}\right)^{(1)} \right) \xi_{\vec\eta}$$ = (42)

$$+ \left(\cal N \rm , \left(\Sigma_{\vec\eta},\Omega_{\vec\eta_\parallel},\Omega_{\vec\eta_\perp}\right)^{(1)} \right)$$

 

$$\dot\xi_{\vec\eta} = i\omega^{(2)}\xi_{\vec\eta} - \left( \frac{\... ...\vec\eta_\parallel},\Omega_{\vec\eta_\perp}\right)^{(2)} \right) \xi_{\vec\eta} - \left( \frac{\partial}{\partial t}{\pmatrix{ \Sigma_{\vec\eta} ... ...vec\eta_\parallel},\Omega_{\vec\eta_\perp}\right)^{(2)} \right) \psi_{\vec\eta}$$ (43)

$$+ \left(\cal N \rm , \left(\Sigma_{\vec\eta},\Omega_{\vec\eta_\parallel},\Omega_{\vec\eta_\perp}\right)^{(2)} \right)$$

 

$$\dot\phi_{\vec\eta} = i\omega^{(3)}\phi_{\vec\eta} - \left( \frac... ...\eta},\Omega_{\vec\eta_\parallel},\Omega_{\vec\eta_\perp}\right)^{(3)} \right),$$ (44)

where we denote $\omega^{(j)} \equiv \omega^{(j)}_{\eta, {\rm max}}$.

Equation (44) is decoupled from the others and its solution is the decrease of the initial value of $\phi$. Equations (42) and (43) are coupled through the linear and nonlinear terms. The interaction through the linear terms is weak, since the coupled linear terms have smaller amplitudes compared to linear terms and they decrease with time due to their dependence on the time derivatives of the eigenvectors, which are very small in the later stages of the shell evolution. The coupling through the nonlinear terms leads to the terms of the third and higher orders, which can be neglected with respect to quadratic terms. Furthermore, the solution of the Eq. (42) has a decreasing character, because the first term on the right side, which includes the "stable'' $\omega^{(1)}$, dominates. Equation (43) is the most interesting one, because it has $\omega ^{(2)}$ in the first linear term, and only the $\omega ^{(2)}$ can be imaginary negative, which has meaning of instability. The explicit form of the Eq. (43) is

 

$$\dot\xi_{\vec\eta} i\omega^{(2)}\xi_{\vec\eta} -\frac{(\vert\dot{\vec\eta}\vert\Sigm... ...rt\vec\eta\vert\Sigma_0}{\vert\vec\eta\vert\Sigma_0(\omega^{(1)}-\omega^{(2)})}$$ =

$$\frac{(i 2\frac{\dot VR-V^2}{R^2} - \dot{\omega}^{(2)})}{\dots} \... ...\vec\eta\vert\Sigma_0\left( \omega^{(1)}-\omega^{(2)}\right) } \sum_{\vec\eta'}$$

$$i\vert\vec\eta-\vec\eta'\vert \Sigma_0 \left( \psi_{\vec\eta-\vec... ...right) \left[ \left( i2\frac{V}{R}-\omega^{(1)}\right) \psi_{\vec\eta'} \right.$$

$$\left. +\left( i2\frac{V}{R}-\omega^{(2)}\right) \xi_{\vec\eta'} ... ...rac{V}{R}} {\vert\vec\eta\vert\Sigma_0\left( \omega^{(1)}-\omega^{(2)}\right) }$$

$$\sum_{\vec\eta'} i \vert\vec\eta'\vert \Sigma_0 \left( \psi_{\vec... ...ft( i2\frac{V}{R}-\omega^{(1)}\right) \psi_{\vec\eta-\vec\eta'} \right. \right.$$

$$\left. \left. +\left( i2\frac{V}{R}-\omega^{(2)}\right) \xi_{\vec... ...] \frac{(\vec\eta-\vec\eta' , \vec\eta')}{\vert\vec\eta-\vec\eta'\vert} \right.$$ (45)

$$\left.+ \phi_{\vec\eta-\vec\eta'} \frac{(\vec\eta-\vec\eta'_{\per... ... +\frac{3V}{R\Sigma_0 \left( \omega^{(1)}-\omega^{(2)}\right)} \sum_{\vec\eta'}$$

$$\vert\vec\eta-\vec\eta'\vert \Sigma_0 \left( \psi_{\vec\eta-\vec\... ...left[ \left( i2\frac{V}{R}-\omega^{(1)}\right) \psi_{\vec\eta'} \right. \right.$$

$$\left. \left. +\left( i2\frac{V}{R}-\omega^{(2)}\right) \xi_{\vec... ..._{\perp} , \vec\eta)}{\vert\vec\eta'_{\perp}\vert\vert\vec\eta\vert} \right\} +$$

$$- \frac{i}{\left( \omega^{(1)}-\omega^{(2)}\right)} \sum_{\vec\et... ...ft( i2\frac{V}{R}-\omega^{(1)}\right) \psi_{\vec\eta-\vec\eta'} \right. \right.$$

$$\left. \left. +\left( i2\frac{V}{R}-\omega^{(2)}\right) \xi_{\vec... ..., \vec\eta')}{\vert\vec\eta-\vec\eta'\vert} + \phi_{\vec\eta-\vec\eta'} \right.$$

$$\left. \frac{(\vec\eta-\vec\eta'_{\perp}, \vec\eta')}{\vert\vec\e... ...a^{(1)}\right) \psi_{\vec\eta'} + \left( i2\frac{V}{R}- \right. \right. \right.$$

$$\left. \left. \left. \omega^{(2)}\right) \xi_{\vec\eta'} \right] ... ...vert\vec\eta\vert} \right\} - \frac{i}{\left( \omega^{(1)}-\omega^{(2)}\right)}$$

$$\sum_{\vec\eta'} \left\{ \left[ \left( i2\frac{V}{R}-\omega^{(1)}... ...{\vec\eta-\vec\eta'} + \left( i2\frac{V}{R}-\omega^{(2)}\right) \right. \right.$$

$$\left. \left. \xi_{\vec\eta-\vec\eta'} \right] \frac{(\vec\eta-\v... ...p}, \vec\eta)}{\vert\vec\eta-\vec\eta'_{\perp}\vert\vert\vec\eta\vert} \right\}$$

$$\left\{ \left[ \left( i2\frac{V}{R}-\omega^{(1)}\right) \psi_{\ve... ...c{V}{R}-\omega^{(2)}\right) \xi_{\vec\eta'} \right] \vert\vec\eta'\vert \right.$$

$$\left. + \phi_{\vec\eta'} \vert\vec\eta'\vert \right\} + \frac{i ... ...(1)}-\omega^{(2)}\right)} \sum_{\vec\eta'} \vert\vec\eta-\vec\eta'\vert\Sigma_0$$

$$\left( \psi_{\vec\eta-\vec\eta'}+\xi_{\vec\eta-\vec\eta'}\right) ... ...+\xi_{\vec\eta'}\right) \frac{(\vec\eta' , \vec\eta)}{\vert\vec\eta'\vert}\cdot$$

In analogy to Fuchs (1996) we group together components with fully imaginary $\omega_{\vec\eta}$ into wave-packets and take the wavenumber $\eta_{\rm max}$. Other components are grouped to the wave-packets of the same width and approximated by the average wave-numbers. The $\xi_{\vec\eta_{\rm max}}$ modes grow, the $\psi_{\vec\eta_{\rm max}}$ and the $\phi_{\vec\eta_{\rm max}}$ modes descend. Other modes oscillate with decrease.

	Figure 3: The evolution of the maximum perturbation of the surface density in the case when the initial values of linear and nonlinear terms of perturbation are in phase.
Open with DEXTER

Demanding $\eta=\eta'=\left\vert\vec\eta-\vec\eta'\right\vert=\eta_{\rm max}$we obtain from geometrical consideration:

$$\begin{array}{ll} (\vec\eta-\vec\eta'_{\pm}, \vec\eta'_{\pm... ...eta'_{\pm}, \vec\eta)=\frac{1}{2}\eta^2_{\rm max} \end{array},$$ (46)

where $\vec\eta'_{+}$ and $\vec\eta'_{-}$ are two wave-vectors inclined at angles $60^{\rm o}$ to the wave-vector $\vec\eta$. It means that $\xi_{\vec\eta}$ wave-packets non-linearly interact with others with wave-vectors $\vec\eta'_{\pm}$. Using Eqs. (29) and (46), the set of Eqs. (45) can be simplified to the form

$$\begin{array}{lclcccr} \dot\xi & = & (i\omega^{(2)} + A) \xi... ...omega^{(2)} + A) \xi_{-} & + & B \xi\xi_{+}^{*} \ , \end{array}$$

where $\xi \equiv \xi_{\vec\eta}$, $\xi_{+} \equiv \xi_{\vec\eta_{+}}$, $\xi_{-} \equiv \xi_{\vec\eta_{-}}$, the asterisk has meaning of the complex conjunction and

$$A =-\frac{-i18RVc^2\dot\Sigma_0\gamma + 18V^2c^2\dot\Sigma_0 +2R^... ...2\Sigma_0\gamma + 9V\dot Vc^2\Sigma_0} {2\Sigma_0(9V^2c^2+\pi^2G^2\rho_0^2R^4)}$$ (47)

and

$$B=\frac{i \pi^3 G^3 \Sigma_0^3 R - i 12 \pi G \Sigma_0 c^2 V \left( \frac{V}{R} + i\gamma\right)} {4c^4\gamma}\ ,$$ (48)

where

$$\gamma=\sqrt{-\frac{V^2}{R^2}-\frac{\pi^2G^2\Sigma_0^2}{c^2}} \cdot$$ (49)

The set of Eqs. (47) describe the time evolution of one triplet of the most interacting modes $(\vec\eta,\vec\eta_{+},\vec\eta_{-})$. A is coming from the time derivative of the eigenvector in (43). Its amplitude is much less than the amplitude of $\omega ^{(2)}$ except around the time $t_{\rm b}$, when $\omega ^{(2)}$ is close to zero. The resulting effect is the time evolution of $\xi, \xi_+, \xi_-$, which for $t > t_{\rm b}$ initially decrease and start to grow only with some delay (see Fig. 4). Nevertheless, the transformation of $\xi, \xi_+, \xi_-$ to $\Sigma$ gives the continuous increase of $\Sigma$ after $t_{\rm b}$ (see Fig. 3). A is the term originating in the transformation of $\Sigma$, $\Omega$ to $\xi, \xi_+, \xi_-$, and it is important around $t_{\rm b}$ only.

	Figure 4: The solution of the set of Eqs. (47) with the initial conditions corresponding to Fig. 3.
Open with DEXTER

6.1 The numerical solution

The set of Eqs. (47) can be solved numerically. We start at the time $t_{\rm b}$ which is the time when the instability begins (imaginary part of $\omega_{\eta,{\rm max}}^{(2)}$ starts to be negative). First we select real and imaginary parts of all initial perturbation amplitudes $\xi, \xi_+, \xi_-$, which have the meaning of initial perturbations of the surface density and of the velocity, such that they correspond to $\Sigma_1 / \Sigma_0 = 0.05$. The magnitude of these perturbations in physical values can be computed from the eigenvectors (31).

The solution is determined by parameters of two types: the first ones, as speed of sound c in the shell, are constant values, the second ones, as the radius of the shell R(t), the expansion velocity V(t) and its time derivative and the surface density $\Sigma_0(t)$ and its time derivative, are functions of time. We can get them either from the analytical Sedov solution (4-6), or from the numerical simulations of the expanding HI shells described by Ehlerová et al. (1997). In this paper we use the Sedov solution (Eqs. 4-6) with following parameters: total energy $E_{\rm tot} = 10^{53}~{\rm erg}$, density of ambient medium $n_0 = 1~{\rm cm}^{-3}$, average molecular weight $\mu = 1.3$, sound speed in the shell $c = 1~{\rm km\,s}^{-1}$.

The time evolution of $\Sigma_0$ and $\Sigma_1$ are presented in Figs. 3, 5, and 7 and the corresponding amplitudes of real and imaginary parts of $\xi, \xi_+, \xi_-$ for the first two cases in Figs. 4 and 6. We can distinguish two situations: the linear and non-linear parts of the perturbation are in phase, so that they support each other, which is seen in Figs. 3 and 4, or they are in anti-phase, so that the nonlinear contribution slows down the linear growth of perturbation, as it is visible in Figs. 5 and 6. Contribution of the non-linear terms depends on the shape of the forming fragments, i.e. on the value of the amplitude functions $\xi, \xi_+, \xi_-$. Figures 3-6 show the extreme cases of that contribution. Intermediate cases, keeping the initial value of the perturbation in surface density at the same level, $\Sigma_1 / \Sigma_0 = 0.05$, are given in Fig. 7. We can also see in Fig. 3 that the maximum contribution of nonlinear terms to the value of the perturbed surface density, at the time when $\Sigma_0 = \Sigma_1$, is $\sim$25% of the linear value.

	Figure 5: The evolution of the maximum perturbation of the surface density in the case when the initial values of linear and nonlinear terms of perturbation are in anti-phase.
Open with DEXTER

All the curves start at the instability time $t_{\rm b}$: functions $\Sigma_{\vec\eta}$, $\Sigma_{\vec\eta_{+}}$ and $\Sigma_{\vec\eta_{-}}$ grow, although the absolute value of the appropriate amplitude functions $\xi$, $\xi_{+}$ and $\xi_{-}$ descend during the short time after the instability begins. It is because $\Sigma_{\vec\eta}$, $\Sigma_{\vec\eta_{+}}$ and $\Sigma_{\vec\eta_{-}}$ are connected to the amplitude functions through the eigenvector (31), whose $\Sigma$ part is always growing with time.

The surface density $\Sigma$ at any point of the tangential plane may be computed at any expansion time after $t_{\rm b}$ using eigenvectors (31) and Eq. (18) written for modes $(\vec\eta,\vec\eta_{+},\vec\eta_{-})$. In Fig. 8 we show the distribution of the surface density $\Sigma$ and in Fig. 9 the velocity field of the surface flows $\vec v$ in the tangential plane for $\eta_{\rm max}$ in the case when the initial perturbations have linear terms in phase with the nonlinear terms corresponding to Figs. 3 and 4 at the time $t = 55~{\rm Myr}$. Figures 10 and 11 give $\Sigma$ and $\vec v$ in the tangential plane at the same time for the case when the linear and nonlinear terms of the initial perturbation are in anti-phase corresponding to Figs. 5 and 6. In the former case, the fragments are well defined and the density peaks are separated one from another with deep depressions in $\Sigma$. In the later case, there are high surface density chains with no distinct peaks and we cannot separate individual fragments.

	Figure 6: The solution of the set of Eqs. (47) with the initial conditions corresponding to Fig. 5.
Open with DEXTER

	Figure 7: The evolution of the maximum perturbation of the surface density for intermediate phase-shifts between initial values of the linear and nonlinear terms of perturbation.
Open with DEXTER

We measure the total mass concentrated in one of the well defined fragments shown in Fig. 8. The mass of this fragment is given in Fig. 12 as a function of time. It decreases because the decrease of its size, which is proportional to $\lambda_{\rm max} = {2 \pi R \over \eta_{\rm max}}$, and increases because the accumulation of the ambient medium and surface flows. After $t_{\rm b}$ the resulting mass of the fragment decreases, since the influence of the size shrinking dominates. This happens when the magnitude of ${{\rm d} \lambda_{\rm max} \over {\rm d}t}$ is larger then the amplitude of the surface flows $v_{\rm max}$ (see Fig. 13). The magnitude of ${{\rm d} \lambda_{\rm max} \over {\rm d}t}$ decreases with time and $v_{\rm max}$ increases, and at $t \sim 53~{\rm Myr}$ they are equal. Since then the inflow dominates and the fragment mass growths.

	Figure 8: The distribution of $\Sigma$ in the tangential plane at the time t = 55 Myr for the perturbation shown in Figs. 3 and 4.
Open with DEXTER

	Figure 9: The velocity vectors of the surface flows $\vec v$ corresponding to Fig. 8.
Open with DEXTER

6.2 The mass spectrum of fragments

After $t_{\rm b}$, when the first mode begins to be gravitationally unstable, more and more modes are unstable and the interval of instability growths. In Fig. 14 we give the values of the fragmentation integral $I_{{\rm f}, \eta}(t)$ as defined in (36) as a function of time. This shows at any time the level of development of a fragment with given $\eta$.

	Figure 10: The same as in Fig. 8 at the same time for the perturbation shown in Figs. 5 and 6.
Open with DEXTER

	Figure 11: The velocity vectors of the surface flows $\vec v$ corresponding to Fig. 10.
Open with DEXTER

	Figure 12: The time evolution of total mass in a well defined fragment.
Open with DEXTER

	Figure 13: $\dot \lambda = {{\rm d}\lambda_{\rm max}\over {\rm d}t}$ and $v_{\rm max}$ as functions of time.
Open with DEXTER

	Figure 14: The time evolution of the fragmentation integral $I_{{\rm f}, \eta }$ as given by Eq. (36).
Open with DEXTER

	Figure 15: The mass spectrum of fragments. The straight line is the power law fit of the decreasing part of the spectrum m-1.4.
Open with DEXTER

Mass $m_{\rm frag}$ of a fragment, which is related to $\lambda = {2 \pi R \over \eta}$, may be defined as

$$m_{\rm frag} = \pi \left({\lambda \over 2}\right)^2 \Sigma_0 = {\pi ^3 R^2 \Sigma_0 \over \eta ^2}\cdot$$ (50)

The formation frequency of fragments corresponding to the wavenumber $\eta$is proportional to $I_{{\rm f}, \eta }$, and the surface of the spherical shell of radius R is able to accommodate $R^2/ \lambda^2$ fragments of wavelength $\lambda$. Therefore, the number ${\rm d}N$ of fragments of mass $m_{\rm frag}$formed out of the spherical shell of radius R in gravitational unstable modes with wavenumber $\eta$ from the interval $(\eta, \eta + {\rm d}\eta)$ is

$${\rm d}N = I_{{\rm f}, \eta} \times {R^2 \over \lambda ^2} {\rm d}\eta.$$ (51)

With this equation, we may derive the mass spectrum of fragments ${\rm d}N/{\rm d}m_{\rm frag}$ at some time. It is shown in Fig. 15 at the time $t_{\rm f}$, when $\Sigma_0 = \Sigma_1$. The masses are between a few times $10^3~{\rm {\it M}_{\odot}}$ and $1.5 \times 10^6~{\rm {\it M}_{\odot}}$ with the highest probability peak at $m_{\rm frag} \simeq 10^4~{\rm {\it M}_{\odot}}$. At that time is the radius of the shell almost $1~{\rm kpc}$ with the expansion velocity $\sim$5 kms^-1, and the total mass collected in the shell is $1.14 \times 10^9~{\rm {\it M}_{\odot}}$.

	Figure 16: Dependence of the fragmentation time on the initial perturbation of the surface density. The vertical line gives the time when the value of $I_{{\rm f}, \eta _{\rm max}} = 1$.
Open with DEXTER

The decreasing part of the mass spectrum can be approximated as a power law ${\rm d}N/{\rm d}m_{\rm frag} \sim m_{\rm frag}^{\alpha }$: the fit of the this part of ${\rm d}N/{\rm d}m_{\rm frag}$ gives $\alpha = - 1.4$, which is close to the observed mass spectrum of GMC in the Milky Way: Combes(1991) gives $\alpha = -1.5$. NANTEN survey of the CO emission of the LMC (Fukui 2001) gives steeper slope of $\alpha = -1.9$, which may be explained in the connection to higher level of random velocities in the LMC compared to the Milky Way resulting in the deficit of high mass clouds.

6.3 The time of fragmentation

The evolution of the maximum perturbation of the surface density can be used to determine the fragmentation time $t_{\rm f}$ of the shell. Because at advanced stages the value of the maximum perturbation rises steeply with the time (see e.g. Fig. 3), we define $t_{\rm f}$ as the time, when maximum perturbation of the surface density is equal to the unperturbed value: $\Sigma_1(t_{\rm f}) = \Sigma_0(t_{\rm f})$. Using the fragmentation integral $I_{{\rm f}, \eta }$ we may compare the development level of different fragments at $t_{\rm f}$ (see Fig. 14). We can say that the most frequent fragments are also the most developed, the more massive form only later.

$t_{\rm f}$ depends on the initial conditions of the set of Eq. (47). They correspond to the initial perturbation of the surface density. We can set them to the value typical for the inhomogeneities in the clumpy interstellar medium ( $10^{19}{-}10^{20}~{\rm cm}^{-2}$), which is at $t_{\rm b}$: $0.01 - 0.2 \times \Sigma_0$.

The dependence of $t_{\rm f}$ on the value of the initial perturbation, $\Sigma_1(t_{\rm b})$, is shown in Fig. 16. Fragments form since $45~{\rm Myr}$, for the largest perturbations, to $70~{\rm Myr}$, for the smallest perturbations. The spread in $t_{\rm f}$ for given $\Sigma_1(t_{\rm b})$ is connected to the different shape of the perturbation as shown in Fig. 7. This time may be compared to the fragmentation time $t_{{\rm f}, \eta }$ obtained for $\eta_{\rm max}$ from the linear analysis defined as a time when $I_{{\rm f}, \eta} = 1$.

7 Conclusions

We evaluate the time evolution of perturbations on the surface of an expanding shell. We complement the linear analysis of the gravitational fragmentation process with the inclusion of nonlinear terms, and we compute the time evolution of fragments after the time when the shell starts to be unstable. Some initial perturbations develop into well separated fragments and we estimate the time evolution of the mass of a fragment, the mass spectrum of fragments, and the spread in their formation time. The computed mass spectrum is close to the observed mass distribution of GMC in the Milky Way, but slightly flatter than the mass spectrum of molecular clouds observed in the LMC. This may be related to higher level of random motions in the LMC compared to the Milky Way, which restricts the formation of late time massive fragments and steepen the resulting mass spectrum. Also interesting is that the more massive fragments form at later times of the shell evolution than the less massive fragments. The formation time depends on the value of the initial perturbation: $t_{\rm f} = 45{-}70~{\rm Myr}$. Large density fluctuations shorten this time and thus in the disturbed ISM with large density fluctuations the fragments form sooner than in quiet and smooth ISM where the density fluctuations are small.

Acknowledgements

We would like to thank Burkhard Fuchs and to anonymous referee for valuable comments. This work was inspired by the paper on the fragmentation of uniformly rotating disks by Fuchs (1996). We are also grateful for an enlighting discussions with B. Fuchs in April 1998 and in March 2000 at Star 2000 conference in Heidelberg. The authors gratefully acknowledge financial support by the Grant Agency of the Academy of Sciences of the Czech Republic under the grant No. A 3003705/1997 and support by the grant project of the Academy of Sciences of the Czech Republic No. K1048102.

References

 
• Bicák, J., & Schmidt, B. G. 1999, ApJ, 521, 708 In the text NASA ADS  
• Bisnovatyi-Kogan, G. S., & Blinnikov, S. I. 1982, Astron. Zh., 59, 876 In the text  
• Bisnovatyi-Kogan, G. S., & Silich, S. A. 1995, Rev. Mod. Phys., 67, 661 In the text NASA ADS  
• Combes, F. 1991, Ann. Rev. Astron. Astrophys., 29, 195 In the text NASA ADS  
• Ehlerová, S., Palous, J., Theis, Ch., & Hensler, G. 1997, A&A, 328, 121 In the text NASA ADS  
• Elmegreen, B. G., Elmegreen, D. M. 1978, ApJ, 220, 1062 In the text  
• Elmegreen, B. G. 1989, ApJ, 340, 786 In the text NASA ADS  
• Elmegreen, B. G. 1994, ApJ, 427, 384 In the text NASA ADS  
• Falle, S. A. E. G. 1981, MNRAS, 195, 1011 In the text NASA ADS  
• Fuchs, B. 1996, MNRAS, 278, 985 In the text NASA ADS  
• Fukui, Y. 2001, Extragalactic Star Clusters, IAU Symp. 207, Pucon, Chile, March 12-18, 2001  
• Lindblad, P. O., Palous, J., Lodén, K., & Lindegreen, L. 1997, The kinematics and nature of Gould's Belt - 30 Myr old star forming region, in: HIPPARCOS Venice '97', ed. B. Battrick, ESA SP-420, 507 In the text  
• Ostriker, J. P., & Mc Kee, C. F. 1988, Rev. Mod. Phys., 60, 1 In the text NASA ADS  
• Redman, M. P., Al-Mostafa, Z. A., Meaburn, J., Bryce, M., & Dyson, J. E. 1999, A&A, 345, 943 In the text NASA ADS  
• Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (Academy Press, New York) In the text  
• Strickland, R., & Blondin, J. M. 1995, ApJ, 449, 727 In the text NASA ADS  
• Vishniac, E. T. 1983, ApJ, 274, 152 In the text NASA ADS  
• Vishniac, E. T. 1994, ApJ, 428, 186 In the text NASA ADS  
• Walter, F., & Brinks, E. 1999, AJ, 118, 273 In the text NASA ADS