Instabilities in the Ionization Zones
Around the First Stars^{1}
Abstract
We consider the evolution of the ionization zone around Population III stars with in protogalaxies with at redshifts , assuming that the dark matter profile is a modified isothermal sphere. We study the conditions for the growth of instabilities in the ionization zones. The RayleighTaylor and thermal instabilities develop efficiently in the ionization zones around 2540 stars, while this efficiency is lower for stars with . For more massive stars (), the flux of ionizing photons is strong enough to considerably reduce the gas density in the ionization zone, and the typical lifetimes of stars ( Myr) are insufficient for the growth of instabilities. The gas in a protogalaxy with with a 200 central star is completely ionized by the end of the star’s lifetime; in the case of a 120 central star, only onethird of the total mass of gas is ionized. Thus, ionizing photons from stars with cannot leave protogalaxies with . If the masses of the central stars are 25 and 40 , the gas in protogalaxies of this mass remains essentially neutral. We discuss the consequences of the evolution of the ionization zones for the propagation of the envelope after the supernova explosions of the stars and the efficiency of enrichment of the intergalactic medium in heavy elements.
1 Introduction
Various instabilities can arise in the ionization zones around massive stars, related to the thermodynamics of the gas behind the ionization front and/or shock, as well as the breakout of the ionization front [13]. These instabilities lead to substantial fluctuations of the density and temperature and increase the velocity dispersion of the gas in the ambient interstellar medium (ISM). Hence, subsequent supernovae will occur in the medium with random flows, probably facilitating a rapid breakup of the supernovae envelopes and increasing the efficiency of mixing of heavy elements in the ISM [4]. These effects are not only important for the ISM [5], but are especially interesting for the problem of the primordial enrichment of the first protogalaxies by Population III stars.
The evolution of ionization zones around the first stars has been studied in some detail [69]. The main focus of these papers is the global dynamics of the propagation of the ionization front and the influence of the ionizing radiation of the first stars on the efficiency of photoevaporation of gas in the first protogalaxies, or in other words, on their destruction [10]. Another problem was the calculation of the number of ionizing photons leaving the protogalaxy, or the escape fraction [11, 12], which is important for studies of the possibility of observing the ionization zones [13] and the dynamics of the ionization of the intergalactic medium [14], and, hence, reionization of the Universe. It is evident that more accurate studies of the evolution of the propagation of the ionization fronts require investigation of the dynamics on small scales and the possible growth of instabilities in ionization zones around the first stars. The first attempt to investigate this problem using threedimensional modeling was carried out in [15], which assumed that the initial density profile in the cloud followed a powerlaw, (steeper distributions were also studied, with the aim of investigating the growth of the RayleighTaylor instability); to some extent, this assumption corresponds to the conditions in the first protogalaxies. However, note that the gas profile also depends on the darkmatter distribution in the protogalaxy, and more realistic studies of the efficiency of the growth of instabilities in the ionization zones taking into account the dynamics of gas cooling in the darkmatter potential are necessary.
Here, we study the growth of instabilities in ionization zones around the first massive stars in spheroidal (nonrotating) dwarf protogalaxies with the darkmatter density profile corresponding to a modified isothermal sphere. We used a cosmological model with a term and cold dark matter (CDM model) with , and assumed a relative abundance of deuterium [16].
2 Model of a protogalaxy and numerical methods
This section presents a brief description of the main model parameters and numerical methods applied (more detail may be found in [17, 19]), as well as the initial conditions used.
2.1 Main Parameters
In the model considered, the protogalaxy consists of gas surrounded by a spherically symmetric halo of dark matter. The darkmatter density profile corresponds to a modified isothermal sphere:
(1) 
where and are the radius of the nucleus and the central density, respectively. We assumed a total mass of the protogalaxy (dark halo and gas) of ; for redshift , this corresponds to perturbations in a CDM model with the parameters derived from the threeyear observations of the CMB radiation by the WMAP satellite [16]. Such a protogalaxy has a virial radius pc (for virial relations, see, for instance, [20]). Simple estimates [21, 22] indicate that protogalaxies at redshift cooled efficiently and produced the conditions necessary for the birth of firstgeneration stars. Because of the low efficiency of the fragmentation of the primary gas, we assume that in the lowmass protogalaxies initially appears only one star (see, e.g., [20]), although under certain conditions, in particular in rotating protogalaxies, there may appear a group of stars [19].
Since the minimum temperature of the primordial gas at large redshifts varies from 40 to 200 K (it is determined by the efficiency of the formation of H and HD molecules [21]), the accretion rates onto protostellar cores were higher than for the modern elemental abundances (). Therefore, first generation stars were much more massive than stars of subsequent generations. Numerical models suggest that their masses varied over a wide range of (see, for instance, [20]). The massive stars formed from the primordial gas were much brighter than stars of subsequent generations containing heavy elements, since the chain requires higher temperatures than the CNO cycle. Of the total mass range, we shall study the ionization zones for 25, 40, 120, and 200 stars. We made this choice because these stars have appreciably different characteristics (see the Table), and such stars ultimately explode as supernovae, but do not collapse into black holes, as is the case with the stars [23].
Mass,  Lifetime, Myr  , s 

25  6.459  
40  3.864  
120  2.521  
200  2.204 
2.2 Numerical Methods
The gas dynamics in the protogalaxy model considered is described using the usual gasdynamical equations in cylindrical coordinates in the approximation of axial symmetry (), which were solved numerically using a finitedifference method with an operator splitting technique [25]. A piecewiseparabolic, thirdorder interpolation scheme was used for the gas transport [26].
To obtain the initial distribution of the gas density in the external gravitational potential of the dark matter , the equilibrium equations in cylindrical coordinates were solved numerically. The gas was assumed to be neutral, to have molecular mass , and to be isothermal, with the temperature equal to the virial value . To solve the equilibrium equations, it is necessary to know the initial velocity of the gas , which we obtain from the relation , where is socalled circular velocity of the gas, which is determined by the gravitational potential of the dark matter only (see the virial relation, for example, in [20]), and is a parameter smaller than unity. The method used to solve the equilibrium equations is described in [27]. The iterations were continued until a gasdensity profile corresponding to the mass of gas inside the virial radius was obtained, . We restricted our consideration to spheroidal (nonrotating) protogalaxies.
2.3 Chemical Kinetics and Radiative Transfer
It was assumed in our model that all chemical agents (atoms, ions, molecules) are passive, i.e., they have a velocity field similar to that of the gas. This makes it possible to consider the gas motion only. The chemical kinetics of the primordial gas takes into account the following main components: H, H, H, He, He, He, H, H, D, D and HD. We computed the electron number density by assuming charge conservation. The helium abundance (by mass) was set to . The rates of chemical reactions for collisional and radiative processes were taken from [28]. The chemicalkinetics equations were solved using the backward differencing scheme [29].
We computed the radiative transfer of the ionizing photons using a grid of radial rays with their origin at the center of the protogalaxy. The number of rays was chosen such that each cell of the twodimensional grid was crossed by at least 10 rays. In contrast to the threedimensional case, it was not necessary to split the rays as the distance to the central source grows for the homogeneous twodimensional grid. The ionizing luminosity was isotropically distributed over all rays at the center of the protogalaxy. In each cell, we summed the number of absorbed photons along each radial ray that crosses the cell, then used this number of photons to compute the rates of photochemical reactions and the corresponding heating rate of the gas. The radiative transfer was computed for three spectral ranges: from the threshold for hydrogen ionization to the threshold for single ionization of helium (13.624.6 eV), from the threshold for single ionization of helium to the threshold for double ionization of helium (24.654.4 eV), and from the threshold for double ionization of helium to infinity. The spectral integration assumed that the spectrum of the central star did not change during the transfer process within each of the three spectral ranges, though the relative weights of these ranges changes during the computations.
The method we used was tested many times, both via comparison with the exact results of simple tests for the velocity of propagation of the ionization front in the spherically symmetric case, and via comparison with other codes for cosmology applications [18].
2.4 Initial Conditions
To improve the resolution at the center of the protogalaxy, we applied a nonuniform computational grid [19]. The grid spacing can be controled using a coefficient A: reducing A increases the resolution in inner regions of the computational grid and reduces the resolution in outer regions. We used and a 900900 grid, which provides a physical resolution better than 0.1 pc in the inner 10 pc, which decreases to 1 pc at pc. The nonuniform spacing is similar along both and .
Our computations assumed that the dark matter in the protogalaxy at redshift was already virialized, and that its configuration can already be described by (1). Later, the gas in the protogalaxy cools and the concentration of gas in the central region grows, becoming cm in our computations. Since physical processes that we did not take into account in our model become important at higher gas concentrations, we assumed that a star is born at the center of the protogalaxy when above number density is attained
3 Numerical computations
Let us consider the evolution of the ionization zones around first stars with masses of 25, 40, 120 and 200 in spheroidal (nonrotating) protogalaxies with total masses of (including dark and barionic matter) that are virialized by redshift . Figure 1 shows the distributions of the gas density and temperature around 25, 40, 120 and 200 stars in the central region of a protogalaxy at , where is the mainsequence lifetime of the star (Table 1, [24]). It is clear that the distribution of the gas in the vicinities of the massive stars depends strongly on the mass of the star. A 200 star ionizes and heats the gas not only in the region kpc, but also outside the virial radius for a protogalaxy, as can be seen in the upper part of Fig. 2. Small fluctuations are visible in the distributions of the density and temperature inside the ionization zone of such a star, and the average temperature of the gas is K. The ionization zone for a lower mass star (120 ) is much smaller, kpc; the fluctuations of the gas density and temperature are somewhat higher than in the gas around a 200 star, but remain small. ”Lobes” of hot (rarefied) and cool (dense) gas are visible in the ionization region, which arise due to fragmentation of the envelope formed by the shock and subsequent compression of these fragments by radiation. The formation of the ”lobes” is related to the shadow of the envelope fragments [3].
The ionization zones for 25 and 40 stars are much smaller than the virial radius of the protogalaxy: 30 and 40 pc, respectively. At the same time, a Dfront probably forms around such stars, in contrast to the weak Rfronts that form around more massive stars (120 and 200 ). Figures 1 and 2 clearly show the ionization front itself and a distinct shock propagating with supersonic velocity relative to the gas around the 25 and 40 stars, and a propagating ionization shock around the 120 and 200 stars. The gas density behind the fronts is fairly high, and cooling effects immediately begin to become important, as is manifest by the enhanced density of the gas envelope formed by the ionization front. Figure 2 shows radial profiles of the gas density and temperature along a ray directed from the center of the protogalaxy at an angle of relative to the vertical axis at times . For the ionization zones around the 120 and 200 stars, the density jump in the envelope created by the Rfront exceeds an order of magnitude at ; later, the envelope smears out (it is disrupted). For the zones around 25 and 40 stars, the profiles in Fig. 2 show that the strong Dfront formed before 0.1 (0.4 and 0.65 Myr, respectively, for 25 and 40 stars) and the distribution of gas behind the front does not show any inhomogeneities. However, at a time of 0.3 , density and temperature fluctuations arise behind the front. These fluctuations may be related to two processes with different natures: thermal instability and RayleighTaylor instability.
The growth of instabilities in the ionization zones depends on the mass of the central star. This is obviously due to the different numbers of ionizing photons emitted by the stars. A 200 star emits 40 times more photons than a 25 star, and the numbers of photons emitted pairs 25 and 40 and 120 and 200 stars differ by factors of two to three (see the table). Therefore, the nature of the instabilities in these pairs is expected to be similar, as is seen in Figs. 1 and 2. On one hand, lowmass stars emit fewer ionizing photons, resulting in a lower velocity for the ionizationfront expansion than in the case of massive stars. Thus, if the zone expands into a region with a powerlaw density profile, at similar ages, the gas density at the front of the ionization zone will be higher for lowmass stars, making the timescales for cooling and the growth of RayleighTaylor instability shorter. On the other hand, the lifetimes of lowmass stars are longer, making the time available for the growth of instabilities longer. Therefore, more developed instability is expected in the vicinities of lower mass stars (2540 ), as is clearly visible in Figs. 1 and 2.
To study variations of the physical parameters along the ionization front, let us find the radius corresponding to the maximum density along a ray directed from the center of the protogalaxy at an angle to the vertical axis. We assume that the circle with radius corresponds to the region of the ionization front, which, in general, is close to the real situation (Fig. 2). Next, we identify the cells of the computational grid that are located closest to the radius ; the polar angle is measured from the vertical axis. As a result, we obtain the axial distribution of the gas along a circle with radius .
Figure 3 presents the distributions of the density (lower) and temperature (upper) along circles with radii pc, and , where , for a 40 star at . The shortwavelength perturbations are modulated by longerwavelength ones. Perturbations with different wavelengths probably arise due to the RayleighTaylor instability, which has an increment , where is the wave vector, is the acceleration of the envelope, and is the ratio of the densities behind and ahead of the front. For the radius pc, a gas velocity km s and the maximum wavelength for which at Myr is pc; this value is close to the scale of the largest perturbations, pc, shown in Fig. 2. Shortwavelength perturbations arose earlier, when the radius of the ionized envelope was smaller. The appearance of the smallest perturbations, with pc, is probably related to thermal instability, since the ionization zone is nonstationary and radiative heating is not completely balanced by cooling. The cooling timescale appears to be fairly short, s, and the size of the perturbation is close to pc, comparable to the size of the computational grid cells. Moreover, RayleighTaylor instability can also be significant on small scales, since its increment grows as .
Let now describe the global characteristics of the ionization zones. Figure 4 shows the evolution of the mass fraction comprised by the gas with relative electron number density exceeding 0.1, in the region inside radii of (solid) and 0.1 (dashed); the mass fraction of gas having a relative electron number density of in these same regions (dotted and dashdotted, respectively); and the mass fraction of gas with temperatures higher than K (dash double dotted). For the 120 and 200 stars (thick solid), the mass of ionized gas increases with time, and exceeds 10% of the mass of the protogalaxy at times and 0.5 , respectively, while the radius of the ionization zone exceeds 50 pc (). At these times, the degree of ionization inside the virial radius is either less than 0.001 or greater then 0.1 (the curves are very close), whereas the temperature of the ionized gas exceeds K. The gas in a protogalaxy with a 200 star is appreciably ionized at 0.8 , while only onethird of the gas is ionized in a protogalaxy with a 120 central star at the end of the star’s lifetime. This provides evidence that ionizing photons from stars with cannot leave protogalaxies. If the mass of the central stars is 25 and 40 (thin curves), the gas in protogalaxies with remains essentially neutral, and the ionization zone does not extend behind 0.1 . There exists a significant fraction of the gas mass with degrees of ionization . It is clear that this gas could be associated with gas that is behind the shock, but has passed through the ionization front. This conclusion is supported by the fact that the mass of gas with K is approximately equal to the mass of gas with (the solid and dashdouble dot curves are very close).
To conclude, we note that the efficiency of the ejection of heavy elements by supernovae is probably much higher in protogalaxies with containing stars with than in such protogalaxies containing stars with , since, first, the ionization fronts from more massive stars substantially reduce the average density of the gas in the protogalaxy and, second, the energy released by the explosions of massive stars is a factor of 310 higher than the corresponding energy for less massive stars [23]. Accordingly, the radiative phase of a supernova remnant will start earlier in protogalaxies with lowermass stars, and the development of RayleighTaylor instability in the envelope is likely, so that a significant fraction of the heavy elements will be locked up in the protogalaxy as a consequence.
4 Conclusions
We have considered the dynamical, thermal, and chemical evolution of gas associated with the formation of the ionization zones around first stars with masses of in protogalaxies with masses of at redshift , and investigated the conditions for the development of instabilities in the ionization zones. We have shown the following.

RayleighTaylor and thermal instabilities develop in the ionization zones, which are especially strong around 2540 stars and less important for stars with masses ; if the stars are more massive (), the flux of ionizing photons over the entire lifetime of a star turns out to be sufficient to create a weak ionization Rfront, behind which the instabilities do not have time to grow.

The gas in a protogalaxy with with a 200 star is completely ionized by the end of the star’s lifetime, while only onethird of the gas is ionized in the case of a 120 star. Thus, ionizing photons from stars with are not able to leave protogalaxies with . If the mass of the central star is 25 or 40 , the gas in a protogalaxy with this mass remains essentially neutral.

After the supernova explosion, heavy elements will be efficiently ejected from protogalaxies with containing stars with masses . The ejection efficiency is negligible for stars of lower massses (2540 ); i.e., after the explosions of 2540 stars, a significant fraction of heavy elements will be locked inside the protogalaxy.
5 Acknowledgements
This work was supported by the Russian Foundation for Basic Research (project code 090200933), the Austrian Scientific Foundation (FWF, project code M 1255N16), the Ministry of Education and Science of the Russian Federation (the Departmental Analytical Targeted Program The Development of the Potential of Higher Education (projects RNP2.1.1.5940, RNP 2.1.1/11879, state contract P685). EOV acknowledges support from the Russian Foundation for Basic Research, project codes 110290701, 110201332) and the Dynasty foundation. Yu.A.S. acknowledges the Russian Foundation for Basic Research (project code 110297124). Computations were performed using computer clusters of the Computing Center of the Southern Federal University and the Center for Multiaccess to Computational Resources of the Southern Federal University. E.O.V. thanks V.N. Datsyuk for useful discussions and advise.
Footnotes
 This paper is published in Astronomy Reports, 2012, Vol. 56, No. 7, pp. 564571.
References
 E.T. Vishniac, Astrophys. J. 274, 152 (1983)
 M.M. Mac Low and M.L. Norman, Astrophys. J. 407, 207 (1993)
 R.J.R. Williams, Mon. Not. Roy. Astron. Soc. 310, 789 (1999)
 M.A. de Avillez, M.M. Mac Low, Astrophys. J. 581, 1047 (2002)
 J. Scalo, and B.G. Elmegreen, Ann. Rev. Astron. Astrophys. 42, 275 (2004)
 D. Whalen, T. Abel, M.L. Norman, Astrophys. J. 610, 14 (2004)
 M.A. Alvarez, V. Bromm, P.R. Shapiro, Astrophys. J. 639, 621 (2006)
 B. W. O’Shea, T. Abel, D. Whalen, M.L. Norman, Astrophys. J. Lett. 628, 5 (2005)
 T. Abel, J.H. Wise, G.L. Bryan, Astrophys. J. Lett. 659, 87 (2007).
 T. Kitayama, N. Yoshida, H. Susa, M. Umemura, Astrophys. J. 613, 631 (2004)
 H. Yajima, M. Umemura, M. Mori, T. Nakamoto, Mon. Not. Roy. Astron. Soc. 398, 715 (2009)
 J.H. Wise, and R. Cen, Astrophys. J. 693, 984 (2009)
 T.H. Greif, J.L. Johnson, R.S. Klessen, V. Bromm, Mon. Not. Roy. Astron. Soc. 399, 639 (2009)
 A. Mesinger, G.L. Bryan, Z. Haiman, Mon. Not. Roy. Astron. Soc. 399, 1650 (2009)
 D. Whalen, and M.L. Norman, Astrophys. J. 673, 664 (2008)
 D.N. Spergel, R. Bean, O. Doré, et. al., Astrophys. J. Suppl. Ser. 170, 377 (2007)
 E.O. Vasiliev, E.I. Vorobyov, Yu.A. Shchekinov, Astron. and Astrophys. 489, 505 (2008)
 I.T. Iliev, B. Ciardi, M.A. Alvarez, A. Maselli, A. Ferrara et al. , Mon. Not. Roy. Astron. Soc. 371, 1057 (2006)
 E.O. Vasiliev, E.I. Vorobyov, Yu.A. Shchekinov, Astron. Reps 54, 890 (2010)
 B. Ciardi and A. Ferrara, Space Sci. Rev. 116, 625 (2005)
 Yu.A. Shchekinov and E.O. Vasiliev, Mon. Not. Roy. Astron. Soc. 368, 454 (2006)
 M. Tegmark, J. Silk, M.J. Rees, A. Blanchard, T. Abel, F. Palla, Astrophys. J. 474, 1 (1997)
 A. Heger and S. Woosley, Astrophys. J. 567, 532 (2002)
 D. Schaerer, Astron. and Astrophys 382, 28 (2002)
 M.L. Norman, J.R. Wilson, R. Barton, Astrophys. J. 239, 968 (1980)
 P. Collela and P.R. Woodward J. Comp. Phys. 54, 174 (1984)
 E.I. Vorobyov, U. Klein, Yu.A. Shchekinov, J. Ott, Astron. and Astrophys. 413, 939 (2004)
 T. Abel, P. Anninos, Yu. Zhang, M.L. Norman, New Astron. 2, 181 (1997)
 E. S. Oran and J. P. Boris, Numerical Simulation of Reactive Flow (Elsevier, New York, 1987; Mir, Moscow, 1990).