---

Smoothed particle hydrodynamics simulations of expanding H II regions

I. Numerical method and applications

T. G. Bisbas¹ - R. Wünsch¹ - A. P. Whitworth¹ - D. A. Hubber^1,2,3

1 - School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff CF24 3AA, UK
2 - Institute for theoretical astrophysics, University of Oslo, Pb 1029 Blindern, 0315 Oslo, Norway
3 - Centre of Mathematics for Applications, University of Oslo, Pb 1053 Blindern, 0316 Oslo, Norway

Received 15 December 2008 / Accepted 10 February 2009

Abstract
Context. Ionizing radiation plays a crucial role in star formation at all epochs. In contemporary star formation, ionization abruptly raises the pressure by more than three orders of magnitude; the temperature increases from 10 K to , and the mean molecular weight decreases by a factor of more than 3. This may result in positive feedback, either by compressing pre-existing clouds and rendering them unstable, or by sweeping up gravitationally unstable shells. It may also result in negative feedback (by dispersing residual dense gas). Ionizing radiation from OB stars is also routinely invoked as a means of injecting kinetic energy into the interstellar medium and as a driver of sequential self-propagating star formation in galaxies.
Aims. We describe a new algorithm for including the dynamical effects of ionizing radiation in SPH simulations, and we present several examples of how the algorithm can be applied to problems in star formation.
Methods. We use the HEALPix software to tessellate the sky and to solve the equation of ionization equilibrium along a ray towards each of the resulting tesserae. We exploit the hierarchical nature of HEALPix to make the algorithm adaptive, so that fine angular resolution is invoked only where it is needed, and the computational cost is kept low.
Results. We present simulations of (i) the spherically symmetric expansion of an H II region inside a uniform-density, non-self-gravitating cloud; (ii) the spherically symmetric expansion of an H II region inside a uniform-density, self-gravitating cloud; (iii) the expansion of an off-centre H II region inside a uniform-density, non-self-gravitating cloud, resulting in rocket acceleration and dispersal of the cloud; and (iv) radiatively driven compression and ablation of a core overrun by an H II region.
Conclusions. The new algorithm provides the means to explore and evaluate the role of ionizing radiation in regulating the efficiency and statistics of star formation.

Key words: methods: numerical - radiative transfer - ISM: HII regions

1 Introduction

The influence of the UV radiation emitted by a massive star (or small group of massive stars) on the ambient interstellar medium is of great interest. As this radiation ionizes the surrounding neutral material, it creates a region of hot ionized gas (Strømgren 1939), known as an H II region. The huge contrast in pressure between the ionized gas in the H II region and the surrounding neutral gas causes the H II region to expand and sweep up the neutral gas into a dense shell (Kahn 1954; Spitzer 1978). This dense shell is likely to be prone to many different instabilities (e.g. Elmegreen 1994; Vishniac 1983), and some of these instabilities may trigger fragmentation and star formation. This mode of star formation is usually referred to as the ``collect and collapse'' mechanism (Elmegreen & Lada 1977; Whitworth et al. 1994a,b). Deharveng et al. (2003, 2005, 2008) and Zavagno et al. (2007) have observed several instances where this mechanism appears to be operating.

The structure of the interstellar medium is observed to be extremely irregular and to contain many clouds. As an H II region expands, it may overrun and compress nearby clouds, causing them to collapse. This mechanism is known as ``radiation driven compression'', and has been simulated by several workers. Some of these simulations (Sandford et al. 1982; Bertoldi 1989; Lefloch & Lazareff 1994; Miao et al. 2006; Henney et al. 2008) concentrate on the morphology of the resulting bright-rimmed clouds. In contrast, Peters et al. (2008) address the possibility that the flow of ionized gas off an irradiated cloud might be an important source of turbulence; and Kessel-Deynet & Burkert (2003), Gritschneder et al. (2009), Dale et al. (2005) and Miao et al. (2009) explore the possibility that the collapse of a bright-rimmed cloud might sometimes lead to triggered star formation. These simulations show that it is possible to reproduce the observed morphologies of bright-rimmed clouds (Lefloch & Lazareff 1995; Lefloch et al. 1997; Morgan et al. 2008). However, the evolution of bright-rimmed clouds and the role of expanding H II regions in triggering star formation are still poorly understood. For example, Dale et al. (2007a,b) argue that the main effect of an expanding H II region may simply be to expose stars that would have formed anyway.

A variety of techniques has been used to follow the propagation of ionizing radiation in hydrodynamic simulations. Kessel-Deynet & Burkert (2000) consider each SPH particle to be at the end of an individual ray, and compute the attenuation of ionizing radiation along this ray in order to determine whether the particle at the end is ionized. Dale et al. (2005, 2007a,b) have refined this method so as to make it less compute-intensive. Alvarez et al. (2006) have invoked the HEALPix algorithm in order to distribute rays from the ionizing source in a regular manner. HEALPix (Górski et al. 2005) is a versatile tree structure which distributes pixels over the celestial sphere, to a user prescribed level of refinement. At each successive level of refinement a pixel from the level above is split into four smaller pixels. Each pixel corresponds to the end of a ray and each ray at a given level represents an equal solid angle on the celestial sphere. Abel & Wandelt (2002) have shown how the HEALPix rays can be split adaptively, to increase the level of refinement where this is needed. This adaptive ray splitting has been used by Abel et al. (2007) and Krumholz et al. (2007).

In this paper we simulate the expansion of an H II region using the Smoothed Particle Hydrodynamics (SPH) code SEREN (Hubber et al., in preparation). SPH (Lucy 1977; Gingold & Monaghan 1977) is a Lagrangian numerical method that describes a fluid using particles. The SEREN code is described in Sect. 3. The transport of ionizing radiation uses the hierarchical ray structure of HEALPix, with rays emanating from the ionizing source into the surrounding medium. We walk along each ray with an adaptive step proportional to the local smoothing length. At each step we calculate the attenuation of ionizing radiation by computing the number of recombinations into excited states of H, i.e. we adopt the on-the-spot approximation (Osterbrock 1974). Rays are split wherever the separation between particles becomes less than the separation between rays. This is a robust criterion, which ensures that a ray is split when it encounters a dense region, whilst a neighbouring ray passing through a more rarefied region is not split. Once the ionization front is located along a ray, the propagation of that ray is terminated, in order to reduce the computational overhead; this is particularly important during the early stages of evolution, when the Strømgren sphere only involves a small fraction of the total number of SPH particles.

Although heavy elements play an important role in determining the temperature of interstellar gas, we do not consider their chemistry here. Instead, we assume that the composition of the gas is X=0.7 hydrogen and Y=0.3 helium, by mass; and that the helium is everywhere neutral. In the neutral gas we assume that the temperature is and the hydrogen is all molecular; hence the mean molecular weight is and the isothermal sound speed is . In the ionized gas we assume that the temperature is and the hydrogen is all ionized; hence and . In the transition region between the molecular and ionized regions, we impose a linear temperature gradient between these two limiting values.

We explore three cases. (i) In the first case, the star lies at the centre of a spherical cloud, and the H II region is spherically symmetric throughout its expansion. (ii) In the second case, the star is placed off-centre inside a spherical cloud; here the H II region breaks out of the cloud on one side, and the remainder of the cloud is accelerated by the rocket effect. (iii) In the third case, the star is located outside the cloud from the outset, and a shock is driven into the cloud ahead of the ionization front (i.e. radiation driven compression). These cases are presented only as illustrative examples of what the code can simulate. Detailed investigations of these phenomena will be presented in subsequent papers.

The paper is organized as follows. In Sect. 2 we discuss briefly the expansion of an H II region, once the initial Strømgren sphere has formed. In Sect. 3 we describe in detail how we treat the propagation of ionizing radiation. In Sect. 4 we test the algorithm on the three cases described above, and in Sect. 5 we discuss the results and conclude.

2 The D-type expansion of an H  II region

2.1 Spitzer solution

Consider an arbitrary density field, , and suppose that there is an ionizing star at the centre of co-ordinates. Assuming ionization equilibrium, and neglecting the diffuse radiation field, the position of the ionization front (IF), in the direction of the unit vector  , is given by , where  is obtained from

Here,   is the mass associated with one hydrogen nucleus when account is taken of the contribution from helium,  is the proton mass, , is the rate at which the exciting star emits Lyman continuum photons, and is the recombination coefficient into excited states only. Equation (1) assumes that the material inside the H II region is fully ionized.

Thus Eq. (1) determines the radius  at which all the ionizing photons emitted in the direction  have been used up balancing recombinations into excited states. We ignore recombinations straight into the ground state by invoking the on-the-spot approximation (Osterbrock 1974).

Strømgren (1939) was the first to show that the transition from a state of almost completely ionized material to a state of almost completely neutral material occurs in a very short distance compared with the dimensions of the H II region. For example, for a spherically-symmetric H II region expanding into a cloud of uniform density  , the radius of the ionization front is given by the Strømgren radius

and the distance over which the degree of ionization changes from to is given by

Here   is the mean photoionization cross section presented by a hydrogen atom to Lyman continuum photons from an OB star.

Because the squared sound speed in the ionized gas inside the H II region is , whereas the squared sound speed in the neutral material outside the H II region is , there is a large pressure difference between the two regimes (more than three orders of magnitude), and this results in rapid expansion of the H II region. The outward propagation of the ionization front is subsonic relative to the ionized gas (where the sound speed is ), but supersonic relative to the neutral gas (where the sound speed is ). Consequently a strong shock front appears ahead the ionization front (for a detailed discussion see Kahn 1954). Spitzer (1978) has obtained an approximate analytic solution for this phase of the evolution. In this solution, the radius of the ionization front is given by

the density of the ionized gas by

the mass of ionized gas by

and the speed at which the ionization front propagates by

2.2 Semi-analytic approach

A similar result can be obtained by solving the equation of motion of the accreting thin shell accelerated by the thermal pressure of the H II region (see e.g. Hosokawa & Inutsuka 2006). We use this approach to obtain a more general solution which includes effects of the gravity and use this solution to calibrate our numerical code with the self-gravity switched on (see Sect. 4.2).

The equation of motion of the shell of mass and radius   is

where , and are the pressure, the mass and the density of the ionized gas, respectively. The second term on the rhs represents the gravitational force due to the H II region mass, and the third term is the average gravitational force acting on the shell due to its own mass (it is zero at its inner edge and at the outer one; see Whitworth & Francis 2002).

The mass of the shell is a difference between the mass of the original neutral gas which occupied the sphere of the radius  and the mass of the H II region

Inserting (9) into (8) and using Eq. (2) we get

We solve this equation numerically using the subroutine integrate.odeint from the Scientific python library (Jones et al. 2001). Since the model assumes an infinitesimally thin shell, the shell radius  can be compared to the average radius of the ionization front and the average radius of the shock front  in the numerical simulations (see Fig. 11 and related discussion in Sect. 4.2).

3 Numerical treatment

We use the SPH code SEREN which has been extensively tested using a wide range of standard and non-standard problems. SEREN is an astrophysical SPH code parallelized using OpenMP and designed primarily for star and planet formation problems. It contains both the traditional SPH method (Monaghan 1992) and the more recent ``grad-h'' SPH formulation (Price & Monaghan 2004). There exists the option to include additional features within the basic SPH algorithm, such as time-dependent and/or Balsara-switched viscosity (Morris & Monaghan 1997; Balsara 1995), artificial conductivity (Price 2008) and a variety of equations of state. The SPH equations can be solved with a choice of integrators, such as the symplectic 2nd order Leapfrog integrator, in conjunction with a block time-stepping scheme.

SEREN uses a Barnes-Hut octal-spatial tree (Barnes & Hut 1986), both to speed up the calculation of gravitational forces, and to facilitate the collation of particle neighbour lists. The Barnes-Hut tree calculates the gravitational force up to quadrupole or octupole order, and has the option of using either the traditional geometrical opening angle criterion or the GADGET-style multipole acceptance criterion (Springel et al. 2001). For problems where multiple gravitationally bound objects form, SEREN uses sink particles (Bate et al. 1995) to follow the simulation beyond the formation of the first dense bound objects. A 4th-order Hermite integrator (Makino & Aarseth 1992) is included, so that the code can follow the later N-body evolution of the sinks, after the accretion phase has finished. A paper describing SEREN in detail, and the tests to which it has been subjected, will be submitted shortly (Hubber et al., in preparation).

3.1 Ray casting

In order to determine the overall shape of the ionization front we consider many different directions  . These directions are chosen using the HEALPix algorithm (Górski et al. 2005). HEALPix generates a set of   directions (hereafter ``rays''), where is the level of refinement and is the lowest level. The rays of level  are distributed uniformly over the celestial sphere, as seen from the ionizing star. Each ray on level  intersects the celestial sphere at the centre of an approximately square element of solid angle

Therefore, the angle between neighbouring rays is of order

When HEALPix goes to a higher level, , the angular resolution is increased by splitting each of the rays at the previous level, , into four rays (see Sect. 3.2). In the simulations presented here, we limit the ray splitting to , so there is a maximum of 196 608 rays at the ionization front.

The integral on the left-hand side of Eq. (1) must be evaluated numerically along each ray until it equals  . To perform this integration we define a set of discrete evaluation points along each ray. At each evaluation point, j, we determine the density  using an SPH summation over the nearest SPH particles, k,

The integral in Eq. (1) is then evaluated using the trapezium method

where f₁ is a dimensionless tolerance parameter of order unity controlling the integration step, and h_j' is the smoothing length of evaluation point j'. The identifier j'=0 denotes the evaluation point placed at the position of the star, and consecutive evaluation points are then located according to

Thus f₁h_j' represents the adaptive step along the line of sight from the ionizing star towards the ionization front. Figure 1 shows a representation of the method. Acceptable accuracy is obtained with f₁=0.25 and this is the value we use in the simulations presented here.

3.2 Ray splitting

The use of ray splitting allows us to maintain good resolution while achieving significant speed-up in comparison with uniform ray tracing (cf. Abel & Wandelt 2002). In our scheme a ray is split into four child-rays as soon as its linear separation from neighbouring rays, , exceeds f₂h_j. Here r_j is the distance from the ionizing star to evaluation point j, f₂ is a dimensionless parameter controlling the angular resolution of the ensemble of rays representing the ionizing radiation, and h_j is the smoothing length of the local evaluation point j. Hence we require

Acceptable results are obtained with f₂ in the range (1.0,1.3). A smaller f₂ gives greater accuracy, but at the expense of more HEALPix rays. As one moves outwards away from the central star, rays are only ever split, they are never merged.

The NESTED version of HEALPix allocates identifiers to rays according to the scheme illustrated in Fig. 2, which shows a small patch on the celestial sphere. Here the mother-ray of ray m has ID

and the four child-rays of ray m have IDs

Here n=0 corresponds to the child-ray to the South, n=1 corresponds to the child-ray to the West, n=2 corresponds to the child-ray to the East and n=3 corresponds to the child-ray to the North. This numbering scheme makes it very simple to trace a ray back to the star at the origin, provided you know the ID of the ray and its level.

3.3 Ray rotation

To avoid numerical artifacts due discretization in the directions of rays, we rotate the ensemble of HEALPix rays through three random angles (about the z-, x'-, and z''-axes), each time we build the ensemble. This process is necessary to prevent the formation of artificial corrugations in the ionization front, as explained in detail in Krumholz et al. (2007).

3.4 Setting temperatures

At each evaluation point we estimate the value of I(r_j) using Eq. (12). If the calculation returns , then the evaluation point lies in the interior of the H II region. Therefore, some of the particles that belong to its neighbour list should be ionized. We flag a particle as ionized, if it passes the following two tests: (a) its coordinates are inside the solid angle of the ray; and (b) its distance from the star is smaller than the distance of the evaluation point (r_j). If a particle is ionized, it is given a temperature . Otherwise it remains neutral with temperature . Figure 3 illustrates the location of the neighbours of evaluation point j, and the ones which pass these tests.

3.5 Location of the ionization front

We store in memory the value of the integral in Eq. (12) and the coordinates of the evaluation point at the end of a given ray segment. The distance  of the ionization front from the star along a given ray is given by the condition . It is determined using a binary chop algorithm between the last two evaluation points r_j-1 and r_j for which . The binary chop algorithm is stopped once the adjustment to r_j becomes smaller than  10^-3 h_j-1.

If the material of the cloud along the line of sight of a ray is sufficiently rarefied, may never reach within the computational domain. We therefore estimate the maximum distance  of all gas particles and if we stop calculating the sum. Any ray that satisfies this condition is characterized as open and the position of the ionization front is not determined along such rays (see Fig. 4).

3.6 Temperature smoothing

It is not feasible to resolve the ionization front in a standard SPH simulation of an evolving H II region (see Appendix A). Consequently we have to smooth the temperature gradient artificially across the ionization front. To demonstrate why this is necessary, we model a uniform density spherical cloud consisting of 10⁶ pre-settled SPH particles. The total mass of the cloud is , and the radius is . At the centre of the cloud there is a single star emitting  ionizing photons per second. There are no gravitational forces. We simulate this system in two ways: in Case A, the temperature gradient is unsmoothed; and in Case B, the temperature gradient is smoothed.

In Case A, the temperatures of the SPH particles are set as described in Sect. 3.4, and therefore there is an abrupt decrease in temperature across the ionization front. The simulation is terminated once the ionization front reaches the edge of the cloud. In Fig. 5 the smoothing lengths of all SPH particles are plotted against their distance from the ionizing source at time . Due to the abrupt temperature and pressure change across the ionization front, a gap is formed between the ionized and neutral regions. This gap prevents the ionization of new material from the neutral region during the expansion, and therefore the mass of the H II region does not increase (see Fig. 7). This in turn alters the evolution of the radius of the H II region, as shown in Fig. 8.

In Case B, the temperature gradient across the ionization front is smoothed. Well inside the H II region, and well outside the ionization front, we set the temperatures of the SPH particles as described in Sect. 3.4. However, in a layer of thickness  across the ionization front we linearly interpolate between the temperatures in the ionized and neutral regimes, i.e.

Here is the smoothing length of the local evaluation point on the ionization front, and  is its radius. In this case, the transition between the two temperature regimes is continuous (see Fig. 6). Crucially, this allows new neutral material to be ionized (see Fig. 7), and the expansion of the H II region now follows the analytical solution (Eq. (4)) rather closely.

Several other authors have encountered the same problem in the past. Kessel-Deynet & Burkert (2000) smooth the ionization fraction over a distance of the order of one local smoothing length at the position where the ionization front is located. Krumholz et al. (2007) considers the same problem and constructs a region in which a similar procedure is followed.

Table 1:   Parameters of the test simulations.

3.7 Computational cost

The computational cost of locating the ionization front scales as , where  is the number of ionized SPH particles. This is demonstrated in Fig. 9, which shows the CPU time, , taken by the ionization subroutine, as a function of , for a sequence of simulations performed with different total numbers of SPH particles, ; these are the same simulations as those used to test convergence in Sect. 4.1.

Abrupt increases in occur when the level of refinement in HEALPix, , increases by one, in response to the resolution criterion described in Sect. 3.2. Apart from this we see that the scaling fits very closely, as might be expected.

In this sequence of simulations, the overall cost of the ionization subroutine, over the entire duration of the simulation, scales as , where is the total number of SPH particles.

4 Applications

In order to demonstrate the capabilities of the new algorithm, we present the results of four different applications. (i) In Sect. 4.1 we revisit the simulation described in Sect. 3.6, in which a spherically symmetric H II region expands inside a spherically symmetric, uniform-density, non-self-gravitating cloud. In Sect. 3.6 this simulation was used to establish the need for temperature smoothing. Here we use it to explore the effect of changing the resolution by varying the total number of SPH particles. We show that with 300 000 SPH particles, the simulations are well converged, both with one another, and with the analytic approximation derived by Spitzer (Eq. (4)). (ii) In Sect. 4.2, we repeat this simulation with a much more massive cloud, and include the effect of self-gravity. The evolution is now well described by the semi-analytic solution derived in Sect. 2.2 (see Eq. (10)). (iii) In Sect. 4.3, we consider the same configuration as in Sect. 4.1, except that the ionizing star is now initially located towards one edge of the cloud, and the self-gravity of the gas is neglected. As a consequence of the intrinsic asymmetry in the initial conditions, the H II region bursts out of the cloud on one side, and the remainder of the cloud on the other side is accelerated away by the rocket mechanism (Oort & Spitzer 1955), and then ionized. (iv) In Sect. 4.4 we consider a much lower-mass, uniform-density, self-gravitating cloud which is overrun by an H II region. The ionization front propagating into the cloud is preceded by a shock front, which compresses the cloud, but does not trigger collapse.

In all simulations, the SPH particles are initially positioned randomly, and then settled to produce a uniform-density ``glass''; the temperature of the ionized gas is set to (except in the boundary layer, where the temperature is smoothed to values between these two extremes, as described in Sect. 3.6); the recombination coefficient into excited states only is taken to be   ; and the volume occupied by the ionizing star is neglected.

The parameters of all the simulations are listed in Table 1. In this table, the first column gives the subsection in which the simulation is discussed, M is the mass of the cloud, R is its radius,  is the temperature of the neutral medium, X and Y are the fractions by mass of hydrogen and helium, respectively, and  are the mean molecular weights of the neutral medium and the ionized medium, respectively,  is the rate of emission of ionizing photons from the exciting star, D is the distance of the exciting star from the centre of the cloud,  is the initial Strømgren radius and is the number of SPH particles used in the simulation.

4.1 Spherically symmetric expansion of an H  II region in a non-self-gravitating uniform-density cloud

In this application, we create a uniform-density spherical cloud having mass , initial radius , and hence initial density   . An ionizing source is placed at the centre of the cloud, and emits ionizing photons at a constant rate . Using Eq. (2), the initial Strømgren radius is . The neutral gas is assumed to be at . There are no gravitational forces.

In order to establish that our code is converged, we evolve this configuration using different numbers of SPH particles:   10⁵, 6  10⁵, 10⁶, and 2  10⁶. We terminate the simulations when the ionization front reaches the edge of the cloud. In Fig. 10 we plot the average radius of the ionization front, , against time for all four simulations. We see that the curves converge for .

In Fig. 11 we plot the radius of the ionization front, , and the radius of the shock front, , against time, for the simulation performed with   10⁶ SPH particles. We determine the radius of the shock front by finding the most distant particle from the source which has radial outward velocity and density . The second condition is necessary because particles near the edge of the cloud move outwards due to the thermal pressure gradient there, long before they are overrun by the shock front. The factor 1.1 is to accommodate numerical noise in the SPH estimate of the density of a particle.

In general, we find that the Spitzer solution (Eq. (4)) predicts the radius of the ionization front well, whereas the semi-analytic solution (Eq. (10)) predicts the position of the shock front. The advantage of the semi-analytic solution is that it can also be used to treat situations in which the self-gravity of the gas is important, as we show in the next application.

4.2 Spherically symmetric expansion of an H  II region in a self-gravitating uniform-density cloud

In this application, we simulate a much larger, uniform-density spherical cloud, having mass M = 1.5  , initial radius , and hence initial density   . An ionizing source is placed at the centre of the cloud, and emits ionizing photons at a constant rate . Using Eq. (2), the initial Strømgren radius is . The neutral gas is assumed to be at . We use 10⁷ SPH particles, and evolve the system for  .

In this simulation self-gravity is taken into account, but with the following two modifications. First, once the radius of the shock front has been determined (as described in Sect. 4.1), we neglect the gravitational acceleration of all the SPH particles outside this radius. This is to prevent infall of the undisturbed neutral gas. Otherwise the outer parts of the cloud are already falling quite rapidly towards the centre by the time the shock front reaches them, and this complicates the interpretation of the dynamics. Secondly, we only take account of the radial component of the gravitational acceleration. We do this to suppress fragmentation of the shell, again in order to keep the dynamics of the shell as simple as possible (i.e. spherically symmetric).

In Fig. 12 we plot the position of the ionization front, the position of the shock front, and the semi-analytic approximation. We see that the position of the shock front agrees very well with the semi-analytic approximation.

4.3 Off-centre expansion of an H  II region

In this application, we consider a uniform-density spherical cloud, having mass , initial radius , and hence initial density   . An ionizing source is placed at distance from the centre of the cloud, and emits ionizing photons at a constant rate ; the initial Strømgren sphere therefore has radius . The neutral gas is assumed to be at . The simulation uses   10⁵ particles. There are no gravitational forces. The simulation is terminated at  .

Figure 13 shows three stages in the evolution of the cloud. At (Fig. 13a), the ionization front breaks through the edge of the cloud on the left-hand (near) side. From this moment on, the ionized gas can stream away freely on this side. However, on the right-hand side the ionization front continues to drive an approximately parabolic shock front into the interior of the cloud. By (Fig. 13b), this shock has passed the centre of the cloud. At this stage the swept-up layer of neutral gas starts to break up into small ``cometary knots''. By (Fig. 13c), the shock front has opened up into a more hyperbolic shape, and has reached the right-hand edge of the cloud. The shocked neutral gas has spawned even more cometary knots. Figure 14 shows a close-up of a few of the knots at . At , about  of the gas remains neutral. All this neutral gas is contained in hundreds of small cometary knots, which are steadily being ablated by the ionizing flux.

We believe, but cannot unambiguously demonstrate, that the formation of these cometary knots is due to a combination of hydrodynamical instabilities. In the early stages, the shell should be prone to the Vishniac instability (Vishniac 1983); and later on, when all the gas has been swept up into the shell and is being accelerated, it should be prone to the Rayleigh-Taylor instability. Therefore we expect the shell to break up into small knots. However, the size of the knots in our simulation is just larger than the resolution limit of the code, i.e. individual knots contain 100 SPH particles, and therefore have masses and diameters . It is not possible to establish the reality of the knots by performing convergence tests; if we increase the total number of SPH particles, , the masses of the knots decrease correspondingly. We note that these knots are very reminiscent of those seen in the Helix Nebula (O'Dell & Handron 1996), and we will attempt to model this source in a future publication (cf. Capriotti 1973; O'Dell & Burkert 1996; Capriotti & Kendall 2006).

4.4 Radiation driven compression

In this application, we consider a uniform-density spherical cloud, having mass , initial radius , and hence initial density   . We place an ionizing source a distance away from the centre of the cloud. The source emits ionizing photons at a constant rate   . Hence the number-flux of ionizing photons incident on the near side of the cloud is   . Since , the rays are parallel and the flux is constant. In this one application we assume X=1 instead of X=0.7, i.e. pure hydrogen with and . The neutral gas is assumed to be at . The sound speeds are therefore and . The simulation uses   10⁵ particles, and self-gravity is taken into account. The simulation is terminated at  .

The above initial conditions are chosen to reproduce - as closely as possible - the initial conditions of the 2D model by Lefloch & Lazareff (1994). We note that from the outset the cloud is over-pressured, and left to its own devices would simply disperse on a timescale of  .

Figure 15 shows the evolution of the cloud. The ionizing flux propagates upwards and rapidly boils off the outer layers on the near side of the cloud. At (Fig. 15a) a shock front starts to compress the remaining neutral gas. At the same time, the north hemisphere starts to expand due to the thermal pressure of the atomic hydrogen. At (Fig. 15b) a dense, prolate, approximately ellipsoidal core forms. By (Fig. 15c), the prolate core has semi-major axis and semi-minor axis ; its mass is . It is therefore thermally sub-critical, and does not collapse to form a star. Instead it is steadily ablated by ionization. By (Fig. 15d), the remnants of the cloud start to develop a cometary tail. By (not shown in Fig. 15), the last vestiges of the cloud are ionized, and they are from the ionizing star.

Figure 16 plots the total mass of neutral gas against time. The undulations seen at , and at are acoustic oscillations, excited as the cloud responds to the increase in external pressure.

Gritschneder et al. (2009; hereafter G09) have also simulated radiatively driven compression. However, the cloud that G09 treat (with mass , radius , and temperature ) is much more massive and much colder than the one we treat here ( , , ), and therefore it is much less resistant to compression and more prone to triggered gravitational collapse. Furthermore, in their first two simulations, G09 use a larger ionizing flux than we do. As a consequence, the clouds in their simulations are more strongly compressed than ours, particularly at the leading edge (i.e. the edge exposed directly to the ionizing flux), and are triggered into gravitational collapse. In contrast, our cloud is more mildly compressed, and evolves towards a centrally condensed configuration but never becomes gravitationally unstable.

5 Discussion and conclusions

We have introduced a new technique for treating the propagation of ionizing radiation in SPH simulations of self-gravitating gas dynamics. The method uses the HEALPix algorithm to tessellate the celestial sphere, and solves the equation of ionization equilibrium along the rays associated with each tessera; rays are split hierarchically to produce greater resolution wherever it is required, i.e. to ensure that the resolution of the radiation transfer matches the resolution of the hydrodynamics, locally. This makes the algorithm very computationally efficient.

The algorithm has been incorporated into the new Smoothed Particle Hydrodynamics code SEREN (Hubber et al., in preparation), and tested against a number of known analytic and semi-analytic problems. These tests are presented here. The code follows the expansion of a spherically symmetric H II region in a non-self-gravitating gas (Spitzer solution); the expansion of a spherically symmetric H II region in a self-gravitating gas (new semi-analytic solution described in Sect. 2.2); the rocket acceleration and subsequent ablation of a massive cloud irradiated by an ionizing star on one side (Oort & Spitzer 1955); and the radiatively driven compression of a pre-existing dense core engulfed by an expanding H II region. The code will be used in future to explore the role of H II regions in triggering and regulating star formation, injecting turbulent energy into the interstellar medium, and eroding molecular clouds.

Acknowledgements

R.W. and A.P.W. acknowledge support by the Human Resources and Mobility Programme of the European Community under contracts MEIF-CT-2006-039802 and MRTN-CT2006-035890. A.P.W. also acknowledges support from the Science and Technology Facilities Council of the UK under contract PP/E000967/1. The computations in this work were carried out on Merlin ARCCA SRIF-3 Cluster. The column density plots presented in this paper were prepared using SPLASH (Price 2007).

Appendix A: Resolution requirements to resolve the ionization front

Consider an H II region having uniform density . The SPH particles have a universal mass  . If each SPH particle has  neighbours, then the diameter, d, of an SPH particle (i.e. the diameter of its smoothing kernel) is given by

Since d is in effect the resolution of the simulation, the ionization front can only be resolved if . Using Eq. (3), this requirement reduces to

Combining Eqs. (6) and (5), we see that the mass of the H II region is given by

The number of SPH particles required to model the H II region is therefore

Evidently this is a prohibitive requirement on , all the more so when one allows that the mass of neutral gas is likely to be even greater than the mass of ionized gas, and so this will require additional SPH particles.

References

 
• Abel, T., & Wandelt, B. D. 2002, MNRAS, 330, L53 In the text  
• Abel, T., Wise, J. H., & Bryan, G. L. 2007, ApJ, 659, L87 In the text  
• Alvarez, M. A., Bromm, V., & Shapiro, P. R. 2006, ApJ, 639, 621 In the text  
• Balsara, D. S. 1995, J. Comput. Phys., 121, 357 In the text  
• Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362 In the text  
• Bertoldi, F. 1989, ApJ, 346, 735 In the text  
• Capriotti, E. R. 1973, ApJ, 179, 495 In the text  
• Capriotti, E. R., & Kendall, A. D. 2006, ApJ, 642, 923 In the text  
• Dale, J. E., Bonnell, I. A., Clarke, C. J., & Bate, M. R. 2005, MNRAS, 358, 291 In the text  
• Dale J. E., Ercolano, B., & Clarke, C. J. 2007a, MNRAS, 382, 1759 In the text  
• Dale, J. E., Clark, P. C., & Bonnell, I. A. 2007b, MNRAS, 377, 535  
• Deharveng, L., Lefloch, B., Zavagno, A., et al. 2003, A&A, 408, L25 In the text  
• Deharveng, L., Zavagno, A., & Caplan, J. 2005, A&A, 433, 565 In the text  
• Deharveng, L., Lefloch, B., Kurtz, S., et al. 2008, A&A, 482, 585 In the text  
• Elmegreen, B. G. 1994, ApJ, 427, 384 In the text  
• Elmegreen, B. G., & Lada, C. J. 1977, ApJ, 214, 725 In the text  
• Gingold, R. A., & Monaghan, J. J. 1977, MNRAS, 181, 375 In the text  
• Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 In the text  
• Gritschneder, M., Naab, T., Burkert, A., et al. 2009, MNRAS, 393, 21 In the text  
• Henney, W. J., Arthur, S. J., De Colle, F., & Mellema, G. 2008, MNRAS, submitted [arXiv:0810.1531] In the text  
• Hosokawa, T., & Inutsuka, S.-I. 2006, ApJ, 646, 240 In the text  
• Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python, http://www.scipy.org/ In the text  
• Kahn, F. D. 1954, Bull. Astron. Inst. Netherlands, 12, 187 In the text  
• Kessel-Deynet, O., & Burkert, A. 2000, MNRAS, 315, 713 In the text  
• Kessel-Deynet, O., & Burkert, A. 2003, MNRAS, 338, 545 In the text  
• Krumholz, M. R., Stone, J. M., & Gardiner, T. A. 2007, ApJ, 671, 518 In the text  
• Lefloch, B., & Lazareff, B. 1994, A&A, 289, 559 In the text  
• Lefloch, B., & Lazareff, B. 1995, A&A, 301, 522 In the text  
• Lefloch, B., Lazareff, B., & Castets, A. 1997, A&A, 324, 249 In the text  
• Lucy, L. B. 1977, AJ, 82, 1013 In the text  
• Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141 In the text  
• Miao, J., White, G. J., Nelson, R., Thompson, M., & Morgan, L. 2006, MNRAS, 369, 143 In the text  
• Miao, J., White, G. J., Thompson, M. A., & Nelson, R. P. 2009, ApJ, 692, 382 In the text  
• Monaghan, J. J. 1992, ARA&A, 30, 543 In the text  
• Morgan, L. K., Thompson, M. A., Urquhart, J. S., & White, G. J. 2008, A&A, 477, 557 In the text  
• Morris, J. P., & Monaghan, J. J. 1997, J. Comp. Phys., 136, 41 In the text  
• O'Dell, C. R., & Handron, K. D. 1996, AJ, 111, 1630 In the text  
• O'Dell, C. R., & Burkert, A. 1997, Planetary Nebulae, ed. H. J. Habing, & H. J. G. L. M. Lamers (Kluwer), IAU Symp., 180, 332  
• Oort, J. H., & Spitzer, L. J. 1955, ApJ, 121, 6 In the text  
• Osterbrock, D. E. 1974, Astrophysics of gaseous nebulae (San Francisco: W. H. Freeman and Co.) In the text  
• Peters, T., Banerjee, R., & Klessen, R. S. 2008, Phys. Scr., T132, 014026 In the text  
• Price, D. J. 2007, PASA, 24, 159 In the text  
• Price, D. J. 2008, J. Comp. Phys., 227, 10040 In the text  
• Price, D. J., & Monaghan, J. J. 2004, MNRAS, 348, 139 In the text  
• Sandford, M. T., II, Whitaker, R. W., & Klein, R. I. 1982, ApJ, 260, 183 In the text  
• Spitzer, L. 1978, Physical Processes in the Interstellar Medium (New York: Wiley-Interscience) In the text  
• Springel, V., Yoshida, N., & White, S. D. M. 2001, New Astron., 6, 79 In the text  
• Strømgren, B. 1939, ApJ, 89, 526 In the text  
• Vishniac, E. T. 1983, ApJ, 274, 152 In the text  
• Whitworth, A. P., & Francis, N. 2002, MNRAS, 329, 641 In the text  
• Whitworth, A. P., Bhattal, A. S., Chapman, S. J., Disney, M. J., & Turner, J. A. 1994a, A&A, 290, 421 In the text  
• Whitworth, A. P., Bhattal, A. S., Chapman, S. J., Disney, M. J., & Turner, J. A. 1994b, MNRAS, 268, 291  
• Zavagno, A., Pomarès, M., Deharveng, L., et al. 2007, A&A, 472, 835 In the text

Footnotes

... HEALPix

http://healpix.jpl.nasa.gov

All Tables

Table 1:   Parameters of the test simulations.

All Figures

---

Copyright ESO 2009